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Preface 



This monograph, which is an extension of Dr. Emadzadeh's doctoral thesis, investi- 
gates the different aspects of utilization of X-ray pulsars for navigation of spacecraft 
in space. In our view, the monograph possesses two unique features. First, it pro- 
vides a solid mathematical formulation for the absolute and relative navigation 
problems based on the use of X-ray pulsar measurements. Second, it presents a 
comprehensive framework for the signal processing techniques needed to obtain the 
navigation solution. 

We had several readers in mind when writing this monograph. One such group is 
the body of university students and researchers who work on new space navigation 
techniques. Using X-ray pulsars for navigation is an interesting field in which there 
are many new challenging problems that need to be addressed. Another target group 
comprises people from the industry. Deep space navigation missions, especially the 
ones directed beyond the solar system, have attracted a lot of attention in recent 
years. Employing new space navigation techniques will definitely play a key role 
in making such missions successful. We hope that the monograph will encourage 
more researchers in the area of space navigation to work on X-ray pulsar-based 
navigation. 

The subject matter requires some familiarity with linear systems, probability, 
and estimation theory. The knowledge is generally assumed to be of advanced un- 
dergraduate and graduate level. It will be more beneficial, if the reader proceeds 
through the chapters in sequence. We first provide some basic background knowl- 
edge on pulsars and a literature review on pulsar-based navigation in Chap. 2. Then, 
we present the navigation problems, and develop the X-ray pulsar signal models in 
Chap. 3. Using these models, we formulate and analyze the pulse delay estimation 
problem in Chap. 3. Different pulse delay estimators are proposed in Chaps. 4, 5, 
and 6. Using the presented estimators, Chap. 7 provides a recursive algorithm to ob- 
tain the navigation solution. Concluding remarks and suggestions for future work 
are given in Chap. 8. 

Finally, we acknowledge Dr. A. Robert Golshan for his valuable comments and 
suggestions, which helped us greatly improve the monograph. 

Los Angeles, CA Amir A. Emadzadeh 

October 2010 Jason L. Speyer 
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Chapter 1 
Introduction 



One of the main requirements for any space mission is to navigate the spacecraft. 
Current space navigation methods highly depend on ground-based operations. 
To achieve more autonomy and also to augment the current available navigation so- 
lutions, we are interested in utilization of celestial-based navigation systems. Such 
systems use signals emitted from celestial sources located at great distances from 
Earth. Because of their special characteristics, X-ray pulsars are potential celestial 
candidates for navigation. 

The purpose of this book is to provide a concise treatment of utilizing a new 
navigation system for space missions, which is based on employing X-ray signals 
emitted from pulsars. Specifically, the book starts with reviewing current navigation 
methods being used for space missions. We provide motivations for adopting new 
space navigation approaches, and we suggest employing X-ray pulsars. We present 
the navigation system structure and provide the mathematical tools required to study 
and analyze different parts of the system. We also develop different signal process- 
ing techniques, which are essential to obtain the pulsar-based navigation solution. 

Chapter 2 begins with an overview of current space navigation systems. It ex- 
plains why utilizing celestial-based navigation systems is an interesting option for 
space missions. It shows why, of all available celestial sources, X-ray pulsars are 
suggested to be employed for space navigation. It also provides a brief treatment of 
pulsars and the history of pulsar-based space navigation. 

In Chap. 3, the structure of the proposed X-ray pulsar-based navigation system is 
explained. Mathematical models describing the X-ray pulsar signals are developed. 
The time of arrival (TOA) of received photons is modeled as a non-homogeneous 
Poisson process (NHPP), and the probability density function of the TOAs is pre- 
sented. Also, an effective algorithm is presented for simulation of the TOAs for any 
given pulsar with a known rate function. Additionally, it is explained how using 
epoch folding, photon intensity function can be retrieved by measurement of the 
TOAs. The noise associated with the procedure is also analyzed. Furthermore, the 
effect of imprecise spacecraft velocity information on epoch folding is studied. 

Chapter 4 focuses on mathematical formulation of the pulse delay estimation 
problem. It employs models of the pulsar photon intensity functions on each detec- 
tor. It also shows how to model the differential time between the spacecraft clocks. 
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2 1 Introduction 

Depending on available spacecraft data, the Cramer-Rao lower bound (CRLB) is 
provided for estimation of the unknown parameters. Some numerical examples are 
also presented. 

In Chap. 5, it is proposed to first recover the pulsar photon intensity functions 
through the epoch folding procedure, and then estimate the pulse delay using the 
recovered photon intensities. Based on epoch folding, two different methods are pro- 
posed for the estimation of the pulse delay: (1) Using the cross correlation function 
between the empirical rate function and the true rate function and (2) Minimizing 
the difference between the empirical rate function and the true one through solving 
a nonlinear least squares problem. Asymptotic behavior of the proposed estimators 
and the effect of spacecraft velocity errors on their performance are also studied. 
Computational complexity of the estimators is investigated as well. Finally, analyt- 
ical results are verified via numerical simulations. 

Based on maximum-likelihood (ML) criterion and direct utilization of the mea- 
sured photon time tags, another pulse delay estimator is proposed in Chap. 6. The 
estimator's asymptotic behavior, and the effect of imprecise spacecraft velocity data 
on its performance is studied. Computational complexity of the proposed estima- 
tor is investigated. And to assess the analytical results, numerical simulations are 
performed. 

In Chap. 7, the dynamics models of inertial measurement unit (IMU) and space- 
craft motion are employed by a Kalman filter to obtain the three-dimensional 
navigation solution. The effect of different navigation system parameters on achiev- 
able estimation accuracies is investigated. In particular, the effects of different 
values of IMU uncertainty, measurement noise variance, the number of pulsars used 
for measurement, their geometric dispersion in the sky map, and imprecise space- 
craft velocity data are considered. 

Finally, concluding remarks and suggestions for future work are provided in 
Chap. 8. 



Chapter 2 

Celestial-Based Navigation: An Overview 



2.1 Introduction 

In this chapter, we present an overview of spacecraft navigation using X-ray pulsars. 
In Sect. 2.2, we present a concise treatment of current navigation methods being 
utilized for space missions. Section 2.3.1 describes why employing celestial-based 
navigation techniques is desirable for space missions. We introduce different types 
of pulsars in Sect. 2.3.2. Section 2.3.3 explains why X-ray pulsars are interesting 
candidates to be used for navigation purposes. A short history of pulsar-based navi- 
gation is given in Sect. 2.3.4. 



2.2 Current Spacecraft Navigation Systems 

Most of space vehicle operations, thus far, have relied widely on Earth-based nav- 
igation methods for absolute position determination [1,2]. Methods such as radar 
range and optical tracking are widely used for this purpose [3]. Although a ground- 
based tracking system has the advantage of not requiring an active hardware on 
the spacecraft itself, it does need extensive ground operations and careful analy- 
sis of the measured data in an electromagnetically noisy background environment. 
Also, as a spacecraft moves further away from Earth, its position estimation error 
increases if a radar-based navigation system is used. To achieve the necessary range 
determination, the radar system needs to know the observation station's position 
on Earth accurately. Another limitation is that such a system requires the knowl- 
edge of positional information of the solar system objects [1]. However, even if 
precise information of the radar station and solar system objects is available, the 
vehicle position estimation can only be accurate to a finite angular accuracy. The 
transmitted radar beam, along with the reflected signal, travels in a cone of uncer- 
tainty. This uncertainty degrades the position knowledge of the vehicle as a linear 
function of distance. Alternatively, many space vehicles, traveling into deep space 
or on interplanetary missions, employ active transmitters for orbit determination 
purposes [1]. The spacecraft receives a ping from an observation station on Earth 
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4 2 Celestial-Based Navigation: An Overview 

and retransmits the signal back to Earth. Then, the radial velocity is measured at 
the receiving station by measuring the Doppler frequency of the transmitted signal. 
Although some improvements are achieved in spacecraft navigation utilizing such 
systems, this method still has errors that increase with distance. Early experiments 
using these tracking systems on the Viking spacecraft showed accuracies to about 
50 km in position estimation error for missions to Mars and positional accuracies on 
the order of hundreds of kilometers at the outer planets [2]. 

Another navigation approach is optical tracking. Spacecraft navigation based on 
optical tracking measurements is performed in a similar way as radar tracking [4]. 
This technique is based on the use of the visible light reflected from a vehicle to 
determine its location. For some optical measurements, a photograph needs to be 
taken and the vehicle's position is calculated after analysis of the photograph and 
comparison to a fixed star background. Therefore, real-time measurements using 
such systems typically are not achieved easily. Furthermore, optical measurements 
are limited by environmental conditions. 

As many missions have concentrated on planetary observations, spacecraft navi- 
gation can be done by taking video images of the planet and comparing them to the 
known planetary parameters such as diameter and position relative to the other ce- 
lestial objects. Throughout this procedure, the position of the spacecraft relative to 
the planet can be determined [5]. This requires the vehicle to be within the vicinity 
of the investigated planet. 

To obtain accurate absolute navigation solutions for deep space missions, a com- 
bination of Earth-based radar ranging and on-vehicle planet imaging is typically 
required. This approach still requires human interaction and interpretation of data. 
Furthermore, as radar-ranging errors increase as the vehicle's distance from Earth 
increases, accurate navigation becomes more complex because of the required finer 
pointing accuracy of ground antennas. Additionally, the vehicles that process video 
images for navigation purposes need to have complicated onboard systems, which 
increase their cost. The imaging process also requires the vehicles to be sufficiently 
close to the planets. Therefore, it is necessary to investigate alternative methods that 
could provide a complete, accurate absolute navigation solution throughout the solar 
system, and perhaps eventually the intergalactic regimes. 

For vehicles operating in space close to the Earth, the current Global Position- 
ing System (GPS) can provide a complete autonomous navigation solution [6]. The 
GPS uses a constellation of between 24 and 32 medium Earth orbit satellites that 
transmit precise microwave signals, enable GPS receivers to determine their loca- 
tion, speed, direction, and time. However, these satellite systems have limited scope 
for the operation of vehicles relatively far from Earth. 

For deep space missions, many spacecraft utilize the Deep Space Network 
(DSN). This system is an international network of antennas that supports inter- 
planetary spacecraft missions, and radio and radar astronomy observations for the 
exploration of the solar system and the universe [7]. The network also supports 
selected Earth- orbiting missions. The DSN currently consists of three deep space 
communications facilities placed ^120° apart around the world: at Goldstone, 
in California's Mojave desert; near Madrid, Spain; and near Canberra, Australia. 
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This strategic placement permits constant observation of spacecraft as the Earth 
rotates and helps to make the DSN the largest and most sensitive scientific telecom- 
munications system in the world. 

Although accurate radial position can be determined using DSN, it requires 
extensive ground operations and scheduling to coordinate the observations. Even 
utilizing interferometry, the angular uncertainty still can increase significantly with 
distance. Position accuracies in the order of l-10km per astronomical unit (AU) of 
distance from Earth are achievable using interferometric measurements of the Very 
Long Baseline Interferometer (VLBI) through the DSN [1]. The VLBI is a type of 
astronomical interferometry used in radio astronomy. It allows observations of an 
object that are made simultaneously by many telescopes to be combined, emulating 
a telescope with a size equal to the maximum separation between the telescopes. 
Data received at each antenna in the array is paired with timing information, usually 
from a local atomic clock, and then stored for later analysis on magnetic tape or hard 
disk. At a later time, the data is correlated with data from other antennas similarly 
recorded to produce the resulting image. The resolution achievable using interfer- 
ometry is proportional to the observing frequency and the distance between the 
antennas farthest apart in the array. The VLBI technique enables this distance to 
be much greater than that possible with conventional interferometry, which requires 
antennas to be physically connected by coaxial cable, waveguide, optical fiber, or 
other types of transmission line. 



2.3 Pulsar-Based Navigation 
2.3.1 Why Celestial-Based Systems? 

Autonomous formation flying of multiple spacecraft is an important technology for 
both deep-space and near-Earth applications [8,9]. One of the main requirements 
of a formation flight is accurate knowledge of the relative position and velocities 
between the vehicles. The spacecraft absolute navigation solution is also needed for 
any space mission. Several researchers have shown that the navigation solution for 
aerial and low-Earth-orbit applications can be obtained by utilizing differential GPS 
(DGPS). DGPS is an enhancement to Global Positioning System that uses a net- 
work of fixed, ground-based reference stations to broadcast the difference between 
the positions indicated by the satellite systems and known fixed positions. These sta- 
tions broadcast the difference between measured satellite pseudoranges and actual 
(internally computed) pseudoranges, with receiver stations correcting their pseudo- 
ranges by the same amount. However, for deep space missions or situations where 
GPS is not available, an alternative approach is needed. Employing Earth-based 
navigation systems, such as DSN, is a possibility. But, as mentioned in Sect. 2.2, 
such systems suffer from low performance in situations where long range naviga- 
tion is required. Furthermore, they are highly based on communicating with Earth 
to analyze their data. 
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Because of the aforementioned problems, the need for higher accuracy, and 
the continuing increase in cost of vehicle operations, spacecraft navigation is 
evolving from Earth-based solutions toward more autonomous methods [10, 11]. 
Autonomous operation means determination of a complete navigation solution by 
the spacecraft to guide itself toward its destination without human interaction or 
assistance. An autonomous navigation system internally computes its own naviga- 
tion and guidance information by using onboard sensors. Any deviation from its 
planned path is detected, reported, and corrected without input from the ground 
mission control. These autonomous operations reduce the dependence of space mis- 
sions on human interaction and communication with Earth. To reduce dependence 
of navigation systems on ground-based operations and achieve more autonomy, 
utilizing celestial-based navigation systems is desirable. Another reason for devel- 
oping such novel navigation methods is to augment current systems by employing 
additional measurements to improve their performance. Celestial-based systems 
use signals emitted from celestial sources located at great distances from Earth. 
Of various celestial sources, X-ray pulsars are interesting candidates for use in both 
absolute and relative navigation systems because of their special characteristics. 
These characteristics are explained in detail in the following. 



2.3.2 Pulsars 

Celestial sources have played significant roles in navigation throughout history, 
although the majority of sources used have been fixed, persistent stars with visi- 
ble radiation. Sources that produce signals with variable intensity are referred to 
as variable celestial sources. There are several classes of variable celestial objects 
emitting signals whose intensities vary from radio signals to gamma-ray over the 
electromagnetic spectrum. Of the different variable source types, ones producing 
a uniquely identifiable signal that is periodic and predictable can be utilized in a 
specific manner for navigation purposes. One particular class of variable celestial 
sources having this property is pulsars. Pulsars are rapidly rotating, highly magne- 
tized neutron stars [12]. As the neutron star spins, charged particles are accelerated 
out along magnetic field lines in the magnetosphere. This acceleration causes the 
particle to emit electromagnetic radiation as a sequence of pulses produced and as 
the magnetic axis (and hence, the radiation beam) crosses the observer's line of sight 
in each rotation (see Fig. 2.1). The repetition period of the pulses is simply the rota- 
tion period of the neutron star. Pulsars are observed in the radio, visible, X-ray, and 
gamma-ray bands of the electromagnetic spectrum [13]. 

Radio pulsars are broadband, stellar pulsating radio sources powered by the ro- 
tation of a neutron star, resulting in a great stability in the pulsar period [14]. Over 
1 ,300 pulsars are known [13], and more are being discovered through new research. 
In some pulsars, irregularities (glitches) in their rotational frequency are observed 
every few years, ranging from the order of 10~ 6 s for the Vela pulsar to only 10~ 8 s 
for the Crab pulsar. Although individually emitted pulses from the pulsars vary over 
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Fig. 2.1 Diagram of a pulsar. Photo courtesy of National Radio Astronomy Observatory (NRAO) 



time, the average pulse shape is stable and characterizes the pulsar. Very precise 
models are established for the mean arrival time of pulsars, whose stability out- 
performs even the most precise artificial time bases. Of all pulsars, the most stable 
ones are the millisecond pulsars. Joseph Taylor and collaborators have demonstrated 
that the timing stability of millisecond pulsars is comparable with terrestrial atomic 
clocks [15]. Millisecond pulsars have been detected in the radio, X-ray, and gamma- 
ray portions of the electromagnetic spectrum. Currently, there are 130 millisecond 
pulsars known in globular clusters. Unfortunately, their signal to noise ratio (SNR) 
is considerably lower than that of longer period pulsars. 

X-ray pulsars are grouped in two different categories according to the source 
of energy that powers the radiation: accretion-powered and rotation-powered 
pulsars [17]. 

1. Accretion-powered pulsars are a class of astronomical objects that are X-ray 
sources displaying strict periodic variations in X-ray intensity. The X-ray pe- 
riods range from as little as a fraction of a second to as much as several minutes. 
An X-ray pulsar consists of a magnetized neutron star in orbit with a normal stel- 
lar companion and is a type of binary star system. The magnetic field strength at 
the surface of the neutron star is typically about 10 12 Gauss, over a trillion times 
stronger than the strength of the magnetic field measured at the surface of the 
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Earth (0.6 Gauss). If the magnetic field and rotation axes of the neutron star are 
misaligned then X-ray pulsations are observed. Accretion pulsars are not stable 
timing sources because their period changes over time. More than 30 accretion- 
powered X-ray pulsars have been discovered with periods from 0.069 to 1,413 s. 
2. Rotation-powered pulsars are rapidly rotating neutron stars whose electromag- 
netic radiation is observed in regularly spaced intervals, or pulses. They differ 
from other types of pulsars in that the source of power for the production of ra- 
diation is the loss of rotational energy. For a long time, the Crab pulsar, the most 
luminous rotation-powered pulsar, had been the only pulsar detected at X-ray 
energies. More than 20 rotation-powered X-ray pulsars have since been detected 
[18, 19]. Figure 2.2 depicts a Chandra X-ray image of the Vela rotation-powered 
pulsar (PSR B083345) [20]. 

Pulsars are the original gamma-ray celestial sources. A few years after the dis- 
covery of radio pulsars by astronomers, the Crab and Vela pulsars were detected at 
the gamma-ray band of the electromagnetic spectrum. Pulsars accelerate particles 
with tremendous energies in their magnetospheres. These particles are ultimately 
responsible for the gamma-ray emission seen from pulsars. The Vela pulsar, which 
spins 1 1 times a second, is the brightest persistent source of gamma rays in the sky. 
Yet gamma rays, the most energetic form of light, are few and far between. Even 
Fermi's Large Area Telescope sees only about one gamma-ray photon from Vela 
every 2min. As opposed to a pulsar's radio beams which only content a few parts 
per million of its total power, gamma-rays represent 10% or more. 



Fig. 2.2 Vela Pulsar (PSR 
B083345) X-ray image taken 
by Chandra X-ray 
observatory (Credit: 
The National Aeronautics 
and Space Administration 
(NASA)/PSU/G. Pavlov 
et al.) 
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By the end of 2004 there were about 1,500 radio pulsars known, but only seven 
had been detected in the gamma-ray band. Pulsars tend to have large magnetic 
fields and to be spinning rapidly. The loss of the pulsar's spin energy eventually 
appears as radiation across the electromagnetic spectrum, including gamma-rays. 
Both observations and models indicate that pulsars eventually lose the ability to 
emit gamma-rays as the pulsar's rotational speed slows down. 

With NASA's Fermi gamma-ray space telescope, astronomers can now have a 
better look at pulsars. In two studies published in the July 2, 2009 edition of Science 
Express, international scientists have analyzed gamma-rays from two dozen pul- 
sars, including 16 discovered by Fermi (see Fig. 2.3). Fermi is the first spacecraft 
able to identify pulsars by their gamma-ray emission alone [21]. The new pulsars 
were discovered as part of a comprehensive search for periodic gamma-ray fluctu- 
ations using 5 months of Fermi Large Area Telescope data and new computational 
techniques. In another part of the study, Fermi team examined gamma-rays from 
eight pulsars, all of which were previously discovered at radio wavelengths. Before 
Fermi launched, it was not clear that pulsars with millisecond periods could emit 
gamma rays. Now it is cleared that they do. It has also become clear that, despite 
their differences, both normal and millisecond pulsars share similar mechanisms for 
emitting gamma-rays. 

NASA's Fermi Gamma-ray Space Telescope is an astrophysics and particle 
physics partnership, developed in collaboration with the U.S. Department of En- 
ergy, along with important contributions from academic institutions and partners in 
France, Germany, Italy, Japan, Sweden, and the U.S. [21]. 
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Fig. 2.3 This gamma-ray all-sky map, which is aligned with the plane of the Milky Way Galaxy, 
shows the pulsar positions, with the 16 new pulsars, detected by Fermi Gamma-ray Space 
Telescope, circled in yellow (eight previously known radio pulsars are in magenta) (Credit: NASA/ 
DOE/Fermi LAT Collaboration) 
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2.3.3 Why Use X-ray Pulsars for Navigation? 

Variable sources emitting signals in the radio band are certainly one potential 
candidate that can be used in a navigation system. However, at the radio frequencies 
that these sources emit, i.e., from 100 MHz to a few GHz, antennas with 20 m in 
diameter or larger are required to detect their signals [22,23]. At these wavelengths, 
"dish" style radio telescopes predominate. The angular resolution of a dish style an- 
tenna is a function of the diameter of the dish in proportion to the wavelength of 
the electromagnetic radiation being observed. This dictates the size of the dish that 
a radio telescope needs to have a useful resolution. 

For most space missions, large antennas highly impact the design and cost of 
the operation [24]. Furthermore, because of neighboring sources that emit in ra- 
dio bands and also low signal intensity of radio pulsars, long integration times are 
needed to obtain a signal with acceptable SNR, suitable for use in a navigation sys- 
tem. Similar limitations exist for the visible variable sources. Additionally, there 
are only five isolated pulsars known to emit in the visible band, and all are faint. 
There are also a only few pulsars discovered which emit in the gamma-ray wave- 
lengths. This is another limitation for utilizing the visible and gamma-ray pulsars in 
a navigation system. 

The disadvantages of the radio and visible sources diminish for sources that emit 
in X-ray band. The main advantage of spacecraft navigation using X-ray sources is 
that small sized detectors can be employed [25]. This provides savings in power and 
mass for spacecraft operations. Another advantage of using X-ray sources is that 
they are widely distributed. The geometric dispersion of pulsars in the sky is impor- 
tant to enhance accuracy of three-dimensional position estimation since the observ- 
ability of the source is an important issue. An important complication that must be 
addressed in utilizing an X-ray source in a navigation system is the timing glitches 
in its rotation rates. Of X-ray pulsars, ones that are bright and have extremely stable 
and predictable rotation rates are suitable candidates for the purpose of navigation. 
These sources are usually older pulsars that have rotation periods on the order of 
several milliseconds. Figure 2.4 provides an image of the Crab Nebula and its pulsar 
(PSR B0531+21), which is the brightest rotation-powered pulsar within the X-ray 
band. The Crab pulsar shows high-flux X-ray emissions with known stable period. 
Hence, it can be considered a suitable candidate for use in navigation. 



2.3.4 History of Pulsar-Based Navigation 

The first pulsar was observed in July 1967 by Bell and Hewish. In 1971, Reichley, 
Downs, and Morris proposed using pulsar signals as a clock for Earth-based systems 
[26]. In 1980, details of methods to determine pulse time of arrivals from pulsar 
signals were provided by Downs and Reichley [27]. In the 1980s and 1990s, it was 
demonstrated that several pulsars matched the quality of atomic clocks [12, 15, 16, 
28]. Because of their stability, pulsars were considered as accurate celestial clocks, 
suitable for navigation. 
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Fig. 2.4 Composite optical/X-ray image of the Crab Nebula pulsar showing surrounding nebular 
gasses stirred by the pulsar's magnetic field and radiation. Photo courtesy of NASA 



In 1974, Downs, a member of the telecommunication division of the Jet 
Propulsion Laboratory (JPL), proposed a spacecraft navigation method based upon 
employing radio signals from a pulsar [29]. Using 27 radio pulsars for navigation 
over an integration time of 24 h, in [29], he showed that an absolute position ac- 
curacy on the order of 150 km was attainable. This introductory paper on pulsar 
navigation serves as the original basis for the work of other researchers in the field. 

During the 1970s, pulsars with X-ray signature were discovered that emit sig- 
nals within the X-ray band of 1-20 keV (2.5el7 -4.8el8 Hz). In 1981, Chester and 
Butman proposed using X-ray pulsars as an option to enhance Earth satellite navi- 
gation [30]. Their research showed that by comparing the arrival times of pulses at a 
spacecraft and at the Earth (via an Earth orbiting satellite), a three-dimensional po- 
sition of the spacecraft can be determined. They reported that a day's worth of data 
from a small onboard X-ray detector yielded a three-dimensional absolute position 
accurate to ~ 150 km. 

In 1988, Wallace studied the issues related to using celestial sources that emit ra- 
dio emission, including pulsars, for navigation applications on Earth [3 1]. He stated 
that the existence of other celestial radio sources obscured weak pulsar signals. As 
expected, radio-based systems require large antennas to detect sources, which make 
them impractical tools for spacecraft. Furthermore, the low signal intensity of ra- 
dio sources requires long integration time to achieve an acceptable SNR. Also, the 
small population of radio pulsars in the optical band of the spectrum severely limits 
an optical pulsar-based navigation system. 
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In 1993, as a part of the NRL-801 experiment for the Advanced Research and 
Global Observation Satellite (ARGOS), Wood proposed a comprehensive approach 
to X-ray navigation covering attitude, position, and time. This study employed X-ray 
sources other than pulsars. As a part of the Naval Research Laboratory (NRL) 
development of this study, Hanson produced a Ph.D. thesis in the field of X-ray 
navigation in 1996 [32]. In his work, he studied attitude determination of spacecraft 
using X-ray pulsars. He used practical data from the HEAO-A1 spacecraft. His ap- 
proach was based on counting the number of received photons, fitting the data to 
preknown curves and minimizing the Chi-squared error. He obtained roll estimates 
with error bias equal to 0.32 deg and standard deviation of 0.030° using a single 
detector and an error bias value equal to 0.012° with 0.0075° of standard deviation. 
He also suggested autonomous timekeeping using X-ray sources by employing a 
phase-locked-loop (PLL). 

In 2004, a research group in Spain revealed a study on the feasibility of an abso- 
lute navigation system based on radio and X-ray pulsars [33]. The group developed 
some models of radio and X-ray pulsar signals, it proposed different algorithms for 
timing estimation of these two categories of celestial sources, studied their perfor- 
mance, and reported the possibility of obtaining absolute position accuracies on the 
order of 10 6 meters [33]. 

In 2005, Sheikh, a member of NRL, produced his Ph.D. thesis in the field of 
X-ray navigation [19, 34, 35]. His work was a part of a research called the X-ray 
Navigation (XNAV), which was directed by the Defense Advanced Research 
Projects Agency (DARPA). He proposed a navigation system based on X-ray 
measurements used by an extended Kalman filter (EKF) for three-dimensional 
position estimation [36,37]. 

Woodfork suggested the use of X-ray pulsars for aiding GPS satellite orbit deter- 
mination in his M.Sc. thesis in 2005 [38]. 

In 2009, Emadzadeh proposed a relative navigation algorithm based on use of 
X-ray pulsar measurements in his Ph.D. thesis [39]. He has studied different aspects 
of the signal processing techniques needed to obtain the X-ray pulsar measurements. 
His dissertation is the main reference of the current book. 



2.4 Summary 

This chapter presents a concise overview of current space navigation methods. 
It provides an introduction on different types of pulsars and their characteristics. 
It suggest using X-ray pulsars to navigate in space for situations where current sys- 
tems are not available or it is desirable to augment the current navigation solutions. 
There are two main reasons that make X-ray pulsars interesting candidates for space 
navigation. One is their stable periodic profile, and the other is that relatively small 
size detectors are needed to detect the X-ray pulsar photons. The latter provides a 
huge advantage for the spacecraft design procedure. A brief history of pulsar-based 
navigation and previous research on this field is also discussed. 



Chapter 3 
Signal Modeling 



3.1 Introduction 

In this chapter, we present an overview of the X-ray pulsar-based navigation system. 
We provide the mathematical tools needed to analyze the proposed system. We also 
characterize the epoch folding procedure. 

Section 3.2 describes the navigation problem and the proposed navigation so- 
lution. In Sect. 3.3, we explain how X-ray pulsar detectors work. Based on the 
presented X-ray detection mechanism, we use some mathematical tools to char- 
acterize the pulsar signals in Sect. 3.4. Section 3.5 offers a clarifying remark on 
pulsar signal modeling. We introduce and mathematically formulate the epoch fold- 
ing procedure in Sect. 3.6. The effect of absolute velocity errors on epoch folding is 
studied in Sect. 3.7. An algorithm is provided in Sect. 3.8 for numerical simulation 
of the X-ray photon TOAs. To verify the analytical results, numerical examples are 
presented in Sect. 3.9. 



3.2 Proposed Navigation System Structure 

First, we consider the relative navigation problem between two spacecraft. At the 
end, we explain how the spacecraft absolute navigation problem can be addressed. 
For relative position estimation, both the space vehicles must lock on the same 
X-ray source, and detect the same signal emitted from it. The spacecraft located far- 
ther away from the pulsar, receives a delayed version of the signal which is observed 
by the spacecraft closer to the source. The distance between the space vehicles is 
proportional to the time delay. The proposed approach is to periodically estimate the 
time delay, and then use these estimates as measurements in a recursive algorithm 
for position estimation. As it can be seen in Fig. 3.1, the relative distance projected 
on the direction of the source, Ad, is related to the time delay, t^, by 

Ad = n- Ax = c?d, (3.1) 

where c is the speed of light, Ax the relative position vector, and n is the normalized 
direction vector pointing to the source. 
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Fig. 3.1 Relative position between two spacecraft observing the same X-ray source 
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Fig. 3.2 Structure of the proposed X-ray pulsar-based navigation system 



The pulsar is far away compared to the relative separation of the spacecraft, 
and the vehicles are assumed to be close enough. Hence, they see the source with 
the same direction vector, and this vector is known. All the measurements are per- 
formed in an inertial frame whose origin is the solar system barycenter (SSB). The 
spacecraft must communicate constantly to share and analyze their data. This com- 
munication channel can be used as an additional information source that can be 
incorporated to enhance the navigation solution. For simplicity, it is assumed that 
the cycle ambiguity problem does not arise. In other words, the signal passes the 
projected distance between the spacecraft, Ad, in less than one pulsar cycle. Rela- 
tivistic effects are negligible. The measurement noise variance is selected based on 
the accuracy of the time delay estimates. On each spacecraft, an inertial measure- 
ment unit (IMU) is to provide the acceleration measurements, which are converted 
to the velocity and position data and are utilized by the Kalman filter. The time de- 
lay estimates are in turn taken in as measurements by the Kalman filter to obtain 
the navigation solution. The precision of relative range is based on the geometric 
distribution of the set of employed pulsars over the sky, in much the same way as 
the geometric distribution of GPS satellites determines the geometric dilution of 
precision. The proposed navigation algorithm can be summarized into the following 
stages (see Fig. 3.2). 
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1. Signal detection using X-ray detectors (Chap. 3). 

2. Pulse time delay estimation (Chaps. 4, 5, and 6). 

3. Obtain the navigation solution utilizing a Kalman filter (Chap. 7). 

Note that if the first spacecraft is at the SSB, the navigation problem reduces to 
the absolute navigation problem of one spacecraft. 



3.3 X-ray Detectors 

X-ray detectors are designed based on measuring the time of arrival (TOA) of 
photons when they hit the detecting material [33]. To obtain the X-ay pulsar mea- 
surements precisely, a low power detection system is required that is capable of 
measuring the photon time of arrivals with sub microsecond accuracy. Detec- 
tors must have a large detection area, low background noise, and be lightweight. 
Furthermore, they need to not require heavy and expensive cooling systems. Previ- 
ous satellite missions have used X-ray detectors with gas-based technology. These 
detectors are prone to leaks in the entrance window [25]. Another disadvantage of 
the detector is that, although made of low density gas, they can be rather heavy be- 
cause of the massive frames they need to contain the gas. Because of these problems, 
an alternative technology is much of interest to make X-ray detectors. 

As part of DARPA XNAV program, NRL has developed a new technology for 
the next generation of X-ray silicon-based detectors that achieve all of the aforemen- 
tioned capabilities [25]. These detectors are made of silicon PIN (p-type, intrinsic, 
and n-type) diodes that are reverse-biased with high voltage so that the whole sili- 
con volume does not have any free charge carrier. Therefore, when an X-ray photon 
hits the silicon, it will create a large number of electrons and holes. These elec- 
trons and holes drift under the electric field to the surface. Then, they are collected 
into preamplifiers and generate signals whose amplitudes are proportional to the 
energy deposited by the X-ray photons. These signals are processed by a shaping 
amplifier and a discriminator to provide the time of arrival of photons. The time 
resolution of the X-ray detector is dictated by the drift time of the electrons across 
the thickness of the detector. The drift time in silicon at room temperature is 20 
ns/mm [25]. The thickness of the detectors determines if they can efficiently detect 
higher energy X-rays. To obtain more accurate solutions for X-ray-based navigation, 
many X-ray photons must be collected. Therefore, the energy range of detectors 
must be maximized, meaning they must be made as thick as possible. The thickest 
silicon detectors currently manufactured are 2 mm thick [25]. More details of the 
silicon-based technology, including required hardware for post-processing of data, 
are explained in [25]. 

Based on this photon detection approach, mathematical models describing the 
X-ray pulsar signals are provided in the next section. 
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3.4 X-ray Pulsar Signal 

Let (to,tf) be the observation interval and r o b s =tf — ?o- Furthermore, let f, be the 
TOA of the (th photon and the set {ti,t2, ■ ■ -Jm}, denoted by {?j}f=i, be a random 
sequence in increasing order, 

to<h<t 2 <---<tM<tf, (3.2) 

where the number of detected photons in (to,tf), M, is also a random variable. 
Let t o = and Nq = then 

N, = max{«, t„ < t} (3.3) 

is called a point process, denoted by {N t , t > 0}. Note that N t is the number of 
detected photons in the interval (0,t). A point process {N t , t > 0} is called a 
non-homogeneous Poisson process (NHPP) with a time-varying rate k(t) > if 
it satisfies the following conditions [40]. 

1 . The probability of detecting one photon in a time interval At, is given by 

P(N t+At - N t = 1) = X(t)At when At -> 0. (3.4) 

2. The probability of detecting more than one photon in At is 

P(N t+ At -N t >2) = when At -> 0. (3.5) 

3. Non-overlapping increments are independent, where 

Nt-N s , t>s (3.6) 

is the increment of the stochastic process {N t ,t > 0}. 

For a fixed t, the number of detected photons in (0,f), N t , is a Poisson random 
variable with parameter J A (| )d| [40], that is 

f/'A(^)d^ expf- fx(^)dA 
P(N, =k) = ^ i ±-^ L (3.7) 



and its mean and variance are 



E[N,} = var[N,} 

= /*A(4)d4 

Jo 
±A(t). (3.8) 
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The number of detected X-ray photons in any fixed time interval (t,s), N t —N s , 
is a Poisson random variable as well with parameter f A (£ )d| [40], i.e., 



P(N,-N s = k) 



A($)d| exp -/ *(4)d$ 



Id 



(3.9) 



The probability density function of the photon TOAs is presented by the follow- 
ing theorem [41]. 

Theorem 3.1. The M- dimensional joint probability density function (pdf) of the 
TOAs set, {/<}£ii, denoted by p ({fj}f=i , Af), is given by 



M 



P ({.,}s„«)=r U m "- 1 



M = 0, 



(3.10) 



where 

A±A(t f )-A(t ) 

is called the integrated rate of the Poisson process. 



(3.11) 



Proof. To calculate p({tt}^,M), consider non-overlapping infinitesimal intervals 
of width Af, symmetric about f,-, i = 1,2, ...,M (see Fig. 3.3). The TOAs pdf, satisfies 
the following, 



At] At] 

ri " Ul-^-,«i + ^- 



'l-Ar,/2 



f M +A( M /2 

tM—AtM/2 



Afjif A?m 

, Tm t I'M z— , W H — y~ 



p({T,}ti,M)dTi---dT M . 



(3.12) 



The expression on the left-hand side of (3.12) is the probability of detecting ex- 
actly one photon in each one of the intervals and none outside them. In other words, 



Afj At] 



At M 



,%€ | tu r— , tM ■ 



At M 

2 



■■P(N tl - Ml/z -N t0 = 0)-P(N tl+&tl/2 -N h - Mr/2 = l)- 



P(N tM+ A tM /2~N tM _ AtM/2 = l)-P(N tf -N tM+AtM/2 = 0). 



(3.13) 
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Fig. 3.3 A sample realization of a Poisson point process 
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Using (3.9), the probability of receiving M = photon in the interval (a,b) is 
given by 

P{N b -N a = 0) = exp (-J\(t)dtj . (3.14) 

Therefore, by letting Af, — > for all i, and substituting (3.4) and (3.14) into the 
right-hand side of (3.13), it can be simplified as 

/ Afi Ati\ ( At M A? M 

Ti £ U i — , t\ + — ,...,%€ Um r- , t M 



2 ' 2 / ' ' V 2 ' 2 

M 

= e- A nA(r«)Af„ (3.15) 

i=l 

where A is given in (3.1 1). 

On the other hand, as A?,- — > for all /, right side of (3.12) equals 

fl+M/2 rtM+^M/2 M 

•••/ P ({Ti}^ 1 ,M)dTi---dTM = p({r,-}^ 1 ,M)nM. (3.16) 

Therefore, since (3.15) and (3.16) are equal, the joint pdf p({^}^i)Af) is 
obtained for M > 1, 

M 

p({t i }tL l ,M)=e- A Y[X(t,). (3.17) 

i=l 



The overall rate function A(f) > 0, represents the aggregate rate of all arriving 
photons from the X-ray pulsar and background, 

A(0=Ab+^A(fet(0) (ph/s), (3.18) 

where h((j>) is the periodic pulsar profile, 0det(O the detected phase, and, Ab and A s 
are the known effective background and source arrival rates, respectively [42]. The 
periodic function h(<j>) is usually defined on the phase interval <j> £ [0, 1) (cycle), and 
then its definition is extended to the entire real line by letting h(<p+n) = h((p), where 
n is an integer. In other words, each period of the pulsar intensity is represented by 
one cycle. Furthermore, the function h(<j>) is non-negative and normalized according 
to Jq /z(0)d0 = 1, and min^ h(<j>) = 0. The shape of the pulsar profile and its period 
are known. Figure 3.4 shows the profile of the Crab pulsar (PSR B0531+21) in the 
X-ray band (1-15 keV), created using multiple observations with the USA experi- 
ment onboard the Advanced Research and Global Observation Satellite (ARGOS) 
vehicle [19]. 
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Phase 

Fig. 3.4 Crab pulsar profile obtained by ARGOS 

The detected rate constants Ab and A s are functions of the detector specifica- 
tions. Specifically, they are proportional to the collective area, A, and the detector 
efficiency, r\ , as 



Ab = b ■ 7] ■ A 

Xs = S-T]'A, 



(3.19) 
(3.20) 



where b is the background rate and s is the source flux. 

Phase at the detector referenced to the beginning of the observation, to, consists 
of an initial phase, 0o € [0, 1), and the accumulated phase 



<foet(t) = <h+ / /(OdT, 
Jt 



(3.21) 



where f (t) is the observed signal frequency, which can be decomposed into two 
different components: X-ray source frequency, / s , and the Doppler frequency shift, 
fd(t). In other words, 

fo(t)=f s +fd(t). (3.22) 

The Doppler frequency f&{t) is due to the detector's velocity, v(t), and is 
given by 

v(0 



fd(t)=fs , 

c 

where c is the speed of light. Consequently, the detected phase equals 



(3.23) 



<j>d e t{t) = <h+Mt-to)+ /d(t)dT. 
Jt 



(3.24) 
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Note that the pulsar is located far deep in space (at infinite distance from the 
detector). Furthermore, the observed phase is calculated referenced to the beginning 
of the observation time, and not the time that photons have left the X-ray source. 
Hence, although h(-) is periodic, the phase ambiguity issue does not arise. 



3.4.1 Constant-Frequency Model 

Assuming that the detector's velocity v(t) — v is a known constant, then 

M0 = to + ('-fo)/o. (3-25) 



where 



./o=(l + ^)/ s . (3.26) 



Hence, from (3.1 8), the rate function is 

Ht;<fo,fo) = fa+kh(<fo + (t-to)f )- (3.27) 

This model is used when the spacecraft moves in a constant radial speed and the 
observed pulse frequency, f , does not change over time or its change is negligible. 
Since h(-) is strictly periodic and its argument is an affine function of time, X{t) 
becomes strictly periodic as well. 

3.4.2 Time-Dependent-Frequency Model 

If v(f) is not constant, f&(t) is a time dependent, and 

<j) d (t)= ff d (r)dx (3.28) 

Jt 

is a nonlinear function of time, which results in a quasi-periodic A(?). This is the 
case when the radial speed of the spacecraft, v(t), changes significantly over the 
observation time. Hence, the observed phase at the detector is 

fcet (*) = *> + ('- fo)/. + to (* ) ( 3 -29) 

and the rate function becomes 

A(f;fo,v(r)) =A<, + a,A(to+(*-<b)/. + &('))• (3.30) 

In (3.30), presence of the nonlinear term, 0</(f ), makes the Poisson rate function 
be quasi-periodic. 
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3.5 Discussion 

Using the TOA pdf presented in (3.10), it can be verified that 

p({T,}^ 1 ,M)dTi---dT M = l, (3.31) 

n 

where Q is the sure event, i.e., the event that any number of photons occur at any 
time instant in the interval [<b,f/]. Note 

[p({r i }f =u M)dr 1 ---dtM=J,f P ({T,}^ 1 ,M)dT 1 ---dT M , (3.32) 

JO m= qJQm 

where Qm is the event of receiving M photons at any M different increasing 
time instants {h,t2, ■ ■ ■ ,?m} in the interval [<0,f/]. The probability of event Qm is 
given by 

p({T,}^ 1 ,M)dT 1 ...dT M 
t 

= p[U =ti e(?o,jy),t2 = f 2 e (h,t f ),...,t M = t M e (t M -x,t f )]. (3.33) 

Note that the sequence {fi,?2, ■ ■ ■ ,?m} in (3.33) is in increasing order. Since there 
exists M ! permutations of tu for a fixed M, the probability that M number of ti's 
occur in no special order, is M\ times that of the sequence occurring in increasing 
order [41]. Hence, 

P[ti = h e (to,t f ),t 2 = t 2 e (t ,tf),, ...,t M = t M e (to,tf)} 

= M\- I p({T,}^ 1 ,M)dT 1 ---dT M . (3.34) 

JQ M 

Using (3.10), the left side of (3.34) can be expressed as 

tf ftf rtf / rtf \ M 

/ •••/>({T;}^ 1 ,M)dT 1 ---dTM = e- A / A(T)dT 
r Jr Jt \Jt j 

= e- A A M (3.35) 

Therefore, from (3.34) and (3.35), for a fixed M 



f A M 

J q p({T / }^ ll M)dTi-dT M = — -e- A . (3.36) 
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Now, using (3.32) and (3.36) 

/^({T,}Ei,M)dT 1 ---dT M -e- A f; 



M=0 



A**_ 



:e- A -e A 



= 1. 



(3.37) 



3.6 Epoch Folding 

A feasible question is how to recover the pulsar rate function employing the mea- 
sured photon TOAs. The procedure of recovering the X-ray pulsar intensity function 
from measured photon time tags is called epoch folding. It is performed as follows: 
all the time tags during the observation time are collected. Then they are folded 
back into a single time interval equal to one pulse period, meaning that their mod- 
ulus with one time period is calculated. Afterward, the period duration is divided 
into some equal-length bins and the number of photons in each bin is counted (see 
Fig. 3.5). Finally, the computed photon counts are normalized and the empirical pul- 
sar profile is derived. Epoch folding is mathematically formulated in the rest of this 
section. 

Assume that the observation time consists of Np pulsar periods. In other words, 
Tabs. ~ N p P, where P is the pulsar period. Let c(?,) be the number of detected pho- 
tons in the (th bin whose center is ?,-. Also, assume that the bin size is 7J, seconds 
and the number of time bins in one period is Ny,, i.e., T/, = P/N\,. Using (3.8), for 
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Fig. 3.5 Epoch folding: All [color-coded] photon TOAs are folded back into the first cycle [to, 
to + P], and it is divided to some equal-length bins. Then, the number of photons in each bin is 
counted and normalized 



3.6 Epoch Folding 23 

infinitesimal 7], values, c(t{) is a Poisson random variable with the following mean 
and variance, 

E[c(t t )]=yat[c(ft)] 

= Kti)T b . (3.38) 

The bin size, Tb, must be chosen small enough such that (3.38) is not violated. 
Now, consider the random variable X (f,), named empirical rate function, which rep- 
resents the normalized number of photons resulted from folding back N p epochs of 
time tags into the ith bin, 

i N P 

Uu)^ WT Y J c j {t i ). (3.39) 

iy P 1 b ;= i 

Then, the relation between A(f,) and A(f,) is given in the following theorem [43]. 

Theorem 3.2. Let A(f,), < 1 t < P, be the true rate function, and X(ti) be the 
empirical rate function, obtained from an epoch folding procedure. Then, 

l{t,) = X{t l )+h{t l ) when r obs ^~ (3.40) 

where n(f,) is referred to as the epoch folding noise. The noise, n(ti), is uncorrelated 
and for long observation times, its mean and variance are given by 

E[n(ti)}=0 (3.41) 

Nu 
var[n(t i )] = -±X(t i ). (3.42) 

-*obs 

Proof Consider the average random variable, C(ti), which obtained from folding 
back N p epochs of time tags into (th bin 

i n p 
£('0 = 7ri<7(fc), (3.43) 

p .7=1 

where j is the epoch's index. Then, from (3.38), as N p —> °° (or equivalently 

Tabs -> °°), 

E[C(t,)}=X(ti)T b (3.44) 

vwr[^N p C(t i )]=X(t i )T b . (3.45) 

Therefore, using (3.39) 

E[l( ti )] = X{t t ) (3.46) 
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Fig. 3.6 Two photon TOAs 
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(3.47) 



As a result of (3.46) and (3.47), the empirical rate function can be modeled by 
(3.40), (3.41), and (3.42). 

To prove that w(f,) is uncorrelated, the first and second moments of N t are calcu- 
lated using (3.8), 



E[N t ] = [ A(r)dr 
Jo 



E[N}] = ( / A(T)dr) + / A(r)dT. 



(3.48) 
(3.49) 



Consider two time instants f, and t^, where t^ > f, (see Fig. 3.6). Because non- 
overlapping increments of N t are independent, 



E [N ti N tk ] = E [N tj (N, k - N ti + N tj )] 

= E[N ti (N tk -N ti )]+E[N*\ 

= E{N tl ]E{N tk -N ll }+E{N? i }. 

Then, using (3.48) and (3.49) 

E[N ti N, k ]= /'A(T)dr( f*A(r)dT+ /'/l(T)dT+l 



(3.50) 



(3.51) 



For infinitesimal values of T b , the number of detected photons in ;th bin, i.e., 
c(tj), can be approximated by 



c(ti)~N t . +Tb/2 -N t ._ Tb/2 . 



(3.52) 



Hence, 



E[c(ti)c(t k )} = E[(N ti+Tb/2 -N ti _ Tb/2 )(N tk+Tb/2 -N tk _ Tb/2 )} 
= E[N ti+Tb/2 N tk+Tb/2 ]-E[N t . +Tb/2 N tk - Tb/2 ] 

-E[N t ,-T b /2N tk+Tb /2}+E[N ti -T b /2N tk -T b /2}- (3.53) 
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Calculating each term in (3.53), using (3.51), it results in 

r'k+T b /2 / fti+T b /2 fti-n/2 

E[c(ti)c(t k )]= A(T)dT / A(T)dT- / A(T)dT 

Jt k -T,,/2 \J0 JO 

= / A(r)dr / A(r)dT (3.54) 

Jtk-T b /2 Jn-Tt/2 

which can be approximated by 

E[cft)c(fc)]«#Afe)A(fc) (3.55) 

for infinitesimal values of 7j. 

As a result, the auto-correlation function of A (/,-) can be calculated from (3.39) 



E[X{ ti )k(t k )] 



1 ^ 2 



N p T t 



%Cj(ti)^Cj(t k ) 

7=1 ,/=l 

I \ 2 N p N p 

1 J^^Eicj^c^)} 



= X{ti)X{t k ). (3.56) 

On the other hand, noting that h(t ) is zero mean, (3.40) leads to 

E[X(ti)X(t k )] = X(ti)X(t k )+E[n(ti)n(t k )]. (3.57) 

Then, from (3.56) and (3.57), the auto correlation function of «(//) is given by 

E[n(ti)n(t k )] = 0. (3.58) 



To perform epoch folding, the observed pulsar frequency, f — l/P, must be con- 
stant and known during the observation time. In other words, the detector velocity, 
v, must be a known constant. Since in practice the velocity is not perfectly known, 
the effect of velocity errors on epoch folding is studied in the next section. 



3.7 Epoch Folding in Presence of Velocity Errors 

As mentioned, it is necessary to investigate how much the epoch folding procedure 
is tolerant to possible errors in the spacecraft velocity data. Let the detector velocity 
error be Av. From (3.26), this error causes the observed frequency to change as 

Av 
A/o=/ s — • (3.59) 

c 
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Since P = l// , A/ makes P change by 



AP=- 



which can be simplified using (3.59), as 



A/o 

f 2 

■to 



3 Signal Modeling 



(3.60) 



./; 



AP=-^Av. 

Hence, the observed pulsar period will be 

P = P + AP. 



(3.61) 



(3.62) 



This pulsar period, which is AP seconds off from the correct one, will be used to 
fold back the photon TOAs into one cycle (see Fig. 3.7). Assuming the number of 
bins is still A^, the new bin size is given by 



Furthermore, the number of epochs to be folded back is 

The error AP, can be translated to the phase domain 

AP 



A(j) 



and using (3.61), it can be simplified by 



P 

foAP 



(3.63) 



(3.64) 



(3.65) 



A<p = -— -Av 

CJo 



— Av. 
c 



(3.66) 



Intensity function 




Fig. 3.7 Epoch folding in the presence of the spacecraft velocity error 
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From (3.38), this bias error in phase changes the statistical properties of cAti) as 
follows. 

E[cj(ti)]=™{cj(ti)} 

= X(tt;h + (j-l)A$)f b , j = l,2,---,N p (3.67) 

In other words, the mean of folded back photon counts from the jth cycle is 
shifted by (/'— 1)A0 cycle. Thus, the empirical rate function equals the follow- 
ing [44]. 

I K 

Mt i ) = WT Y jC j{t i ) 

N p T h j=l 



I n p 



N P.i=i 



£A(fc;fo + (j-l)A*)+ftfc). (3-68) 



where 



and 



E[n(t t )]=0 (3.69) 



var[«(f,)] = * J A.(tt',4o + U- W)- ( 3 - 70 ) 

™p l b j=\ 

From (3.70) it is clear that if (j — 1)A0 is always smaller than the bin size, (3.70) 
boils down to the result given in (3.40). That is, if 

|A0|(* p -l)«|A0|# p <i- (3.71) 

or, using (3.66) 

Av<— C — (3.72) 

then, the empirical function satisfactorily represents the true rate function. The up- 
per bound in (3.72) shows how much error on the spacecraft velocity data can be 
tolerated by the epoch folding procedure. If the velocity error exceeds the bound, as 
(3.70) shows, the empirical rate function will be a deteriorated version of the true 
function. The deterioration shows as a shift in the initial phase, and a change in the 
magnitude of the rate function. 

An interesting point (3.72) shows is that the upper bound is inversely proportional 
to the number of bins N\,. But one should be careful in choosing N\, too small, 
because (3.38) must not be violated. 
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3.8 Generating Photon TOAs 

To assess the analytical results, the X-ray photon TOAs corresponding to a particular 
known pulsar rate function X (t) must be simulated. To realize corresponding TOAs, 
an algorithm based on inversion of the integrated rate function is used [43,45]. 
Let U be a uniform random variable in the interval (0, 1). Define 

g( M )=F- 1 («), (3.73) 

where F^T 1 (•) is the inverse function of F y (-). Then, Y — g(U) is a random variable 
with the distribution F y (y) [46]. That is, 

If y = F-\U) then P{Y < y} = F y (y). (3.74) 

Now let F z \t n =t (z\t„ — t) be the probability distribution of Z = t n+ \ — t n given 
t n = t, where ti is the TOA of the ith photon. Then, 

P{Z > z\t„ =0 = 1- F z \ tn=t {z\t n = t). (3.75) 

Also, using (3.9) 

P{Z>z\t n =t)=P{N t+z -N,=0) 

= exp( - / X(t)dt 



,(-(A(t+z)-A(t))\ 



= expl-(A(r + z)-A(f)) I, (3.76) 

where A (t ) is defined in (3.8). Therefore, from (3.75) and (3.76) 

Fz\t„=t(z\tn=t) = l-ex V {-(A(t + z)-A(t))} (3.77) 

and 

K\tl=M\ t n = t ) = - t + A ^{Mt)-W-z),) (3.78) 

where FT 1 (•) is the inverse function of F : u =t (-). Hence, given t n = t , Z can be 
generated as 

Z=-f+A- 1 (A(f)-ln(l -£/)). (3.79) 

As 1 — U is also a uniform random variable in the interval (0, 1), without loss of 
generality, (3.79) can be restated as 

Z = -t+A- l (A{t)~lnU). (3.80) 
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Also, note that if X is an exponential random variable, described by the following 
pdf, 

pWs= /v-4* *>o (381) 

I otherwise 
then, its probability distribution is given by 

F x (x) = l-e~^ eX (3.82) 



and 



Therefore, from (3.74) 



F-\x) = -^-\n{\-x). (3.83) 

A e 



X = -^\n{\-U) (3.84) 

A e 

is an exponential random variable with parameter X e . As 1 — t/ is also a uniform 
random variable in the interval (0, 1), it can be concluded that 

X = - — \nU (3.85) 

A e 

is an exponential random variable with parameter A e . Thus, the term — In U in (3.80) 
is an exponential random variable, named E, with parameter A e = 1. As a result, 
using (3.80), given t — t n , t n+ \ can be generated as 

tn+i =t n +z 

= A- 1 (A{t„)+E). (3.86) 

The preceding discussion is concluded into Algorithm 3.1. 

Algorithm 3.1 TOA Simulation: 

• L <— (L is an auxiliary variable) 

• k <— (k is a counter) 

• WHILE r* < f/ 

- Generate an exponential random variable, E, with parameter \ e = 1. 

- fc<-Jfc+l 

- L«-A- 1 (A(L) + £) 

- t k <-L 

• END 
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3.9 Numerical Examples 

The Crab pulsar (PSR B0531+21) is selected to verify the analytical results. The 
pulsar period is 33.5 ms, and its profile, h(<j>), is shown in Fig. 3.8. The constant 
arrival rates are chosen to be X\, = 5 and A s = 15 ph/s. The detector's velocity is 
assumed to be v = 3 km/s; and the initial phase observed at the first detector is 
00 = 0.16 cycle. 

A Monte Carlo simulation was performed in which 500 independent realization 
of photon TOAs were processed. Photon TOAs were generated as realizations of a 
non-homogeneous Poisson process, as explained in Algorithm 3.1. The accumulated 
rate function used to generate the TOAs, A, is plotted in Fig. 3.8. 

An observation time of T Q ^ = 100 s is selected to verify the epoch folding 
results. The photon time tags are folded back into one single pulsar cycle which 
is divided into Nt, = 1,024 bins, and the rate functions on each detector are de- 
rived using (3.39). The empirical rate function along with the true one is plotted 
in Fig. 3.9. As the plots show, the empirical rate function converges to the true 
function. 

In Fig. 3.9, the variance of h(t), obtained via simulation, is also plotted and 
compared to the analytical variance for each detector. The plots show an excellent 
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Fig. 3.8 Crab pulsar profile, and A (t ) for Ab = 5 and A s = 15 ph/s 
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Fig. 3.9 Pulsar rate function, and the variance of h(t) for r o t, s = 100 s 
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agreement between the simulation and the analytical value, given in (3.42). The 
difference between these two plots gets larger at the time instants where there is a 
sharp change in the rate function. This phenomenon is due to the limited number 
of time bins and finite observation time which violate (3.38). In other words, if the 
observation time is infinitely long and the time bins are infinitely small then these 
jumps in empirical graphs will vanish. 

The auto-correlation of X (t), E[X(t)X(t + Dq)], is calculated numerically for the 
time delay Dq — P/A second, where P is the pulsar observed period. The plot is 
shown in Fig.3.10, and it verifies that E[X(t)X(t +D )] = X(t)X(t+D ), This 
means that the epoch folding noise, n(t), is uncorrected. 

The effect of spacecraft imprecise absolute velocity data on epoch folding is 
studied as well. The velocity error is intentionally chosen to be big: Av = 1 ,000 m/s. 
For the selected parameters, the upper bound (3.72) roughly equals 10 m/s. Since 
the velocity error is bigger than the bound, the empirical rate function and the noise 
variance are expected to be deteriorated as opposed to the case where the veloc- 
ity is perfectly known. Figure 3.1 1 verifies this fact. Comparing the empirical rate 
function to the true one shows that it is shifted, and its peak is shortened and flat- 
tened. The figure also verifies that the variance of the epoch folding noise is shifted 
because of the velocity error. 
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■E[X(t)X(t + P/4)} 
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Fig. 3.10 Autocorrelation of X (t ) for D = P/4 
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Fig. 3.11 X (t) and the variance of h{t) when Av = 1 ,000 m/s 
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3.10 Summary 

In this chapter, we offer an overview of the proposed navigation system structure. 
We explain how by estimation of the time delay between the X-ray pulsar signals, 
the spacecraft relative or absolute position can be estimated. More detailed math- 
ematical formulation of the estimation problem, and our proposed solutions to the 
problem will be presented in the following chapters. We also suggest using IMUs to 
provide the spacecraft acceleration data. We explain how based on time tagging the 
photon time of arrivals, X-ray pulsar detectors work. Using the presented detection 
mechanism, we develop the mathematical models needed to characterize the pulsar 
signals. We formulate and analyze the epoch folding procedure and show how it 
results in recovering the pulsar intensity function. We investigate how errors on the 
spacecraft velocity data can affect the epoch folding, and show that they result in 
recovering a deteriorated rate function. Finally, we examine the theoretical results 
via numerical simulations. 



Chapter 4 

Pulse Delay Estimation 



4.1 Introduction 

As explained in Chap. 3, the navigation system measurement is obtained through 
estimation of the time delay between the received signals. The delay estimation 
problem plays the most important role in the navigation system. In this chapter, we 
study this problem in more detail. 

Section 4.2 offers a mathematical formulation of the delay estimation problem. 
In Sects. 4.3 and 4.4, we present a lower bound on the variance of the estimation 
error, called the Cramer-Rao lower bound (CRLB). Some numerical examples on 
calculation of the CRLB is also given in Sect. 4.5. 



4.2 Pulse Delay Estimation 

Let {t; }:Ji be the measured photon TOAs at the first detector. This TOA set cor- 
responds to the following rate function. 

Xi(t,to)=h + A. h (ft + {t- f )/i ) (4. 1) 

At the second spacecraft, a different TOA set is measured, {t\ } t 2v which cor- 
responds to the delayed rate function observed by the first detector 

h(t;<i)i)=h+Kh((h + {t-to-t e -t x )f2), (4.2) 

where t x is the true time delay between the detected rate functions, and t e is the 
differential time between clocks [47]. Let t& = t e + t x , and define fc as 

fc = fr-/2k. (4.3) 
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Then, the rate function Az(-) in (4.2) can be simplified as 

Mr <h) = k> + ^ h (<h + ( f - f o)/2) ■ (4.4) 

Note that f\ and fc are observed pulsar frequencies. To construct the measure- 
ments for the Kalman filter, t^ must be estimated. 

To estimate the pulse delay, t^, <j>i and (fa will be estimated on each detector. Then, 
using (4.3), t& can be found. We have proposed two different methods for estimation 
of the pulsar intensity function's initial phase, and analyzed their characteristics. 

1 . Recovering the pulsar rate functions at each detector through the epoch folding 
procedure represents the first strategy. The empirical rate function will be used 
for estimation of the initial phase values. Two different pulse delay estimators 
based on epoch folding are introduced, and their performance against the CRLB 
is analyzed. One uses the cross-correlation function between the empirical and 
true photon intensity functions [44]. The other employs a least-squares criterion 
[43,44,48]. Details are given in Chap. 5. 

2. The second strategy is based on the direct use of the measured photon TOAs. To 
estimate the unknown parameters (initial phase values and Doppler frequencies), 
a maximum likelihood (ML) estimation problem is formulated using the pdf of 
the photon time tags [43,49]. Different aspects of this approach are explained in 
Chap. 6. 



4.3 CRLB 

Being able to place a lower bound on the variance of any unbiased estimator is 
extremely useful in practice. It allows one to assert if the estimator is the minimum 
variance estimator and provides a benchmark against which the performance of any 
unbiased estimator can be compared. An estimator that is unbiased and attains the 
CRLB is said to be efficient, meaning that it efficiently uses the data. The CRLB is 
stated in the following theorem [50]. 

Theorem 4.1. Let p(x;6) be the probability density function (or probability mass 
function) of the observed random vector, x, which is parameterized by the unknown 
parameter 9. A semicolon is used to denote the dependence. Assume that p(x;9) 
satisfies the following "regularity" conditions 



^np(x;0) 



0, for allO, (4.5) 



where 6 is the unknown parameter vector, and the expectation is taken with respect 
to p(x;6). Then, the covariance matrix of any unbiased estimator 6 satisfies 



cov(9)>r l {G), (4.6) 
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where the inequality ">" means (cov(Q) — I~ l (6)) is positive semidefinite and 
1(d) is the Fisher matrix given by 



W)]ij = -E 



d 



2 



dOidOj 



\np(x;8) 



(4.7) 



where the derivatives are evaluated at the true value of and the expectation is 
taken with respect to p(x;0). The matrix I (6) is called the Cramer-Rao lower 
bound. 

Now consider the following pulsar rate function, which corresponds to the mea- 
sured TOA set on the detector. 

a.(f;to,/ ) =h + Kh(<fo + (t-t )f o ). (4.8) 

Assuming that 0o and f are not known, the aim is to find the CRLB for es- 
timation of these parameters. The result is presented in the following theorem 
[51,52]. 

Theorem 4.2. Let 

= (fo /o) T (4.9) 

be the unknown parameters vector. The Fisher matrix and its inverse (CRLB) for 
estimation of 9 in (4.8) are given by 

m = k( 6 Y T t), (4.10) 

6 V3r o 2 bs 27i; 



where 



Therefore, 



,W'jyML df ,4,,) 

Jo A b + A s /!(>) 



CRLB(0)=/- 1 (0) 



2 / 2/T ohs -3/r 2 b ; 



Ip\-3/Tl s 6/r c 



(4.12) 



obs 



Proof Let 6>i and 62 be the elements of the vector (this notation is used for 
simplicity). To derive the CRLB, the photon TOA measurements are transformed 
into another set of data containing the number of detected photons. This new 
data set is formed by dividing the observation time into N number of bins whose 
length is At and by placing the number of counted photons in each bin into the 
vector 

T 

x= (x\ X2 ■■■ x N ) (4.13) 
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see* 



*-«k 



■x- 



v N-[ -^N 



x*e<-x-x 



A? 



Fig. 4.1 Binning the observation time 



as shown in Fig. 4.1. For sufficiently small At, the rate of arrival in nth bin can be 
assumed to be constant: 



1 H +nAt 

A„(0)«-/ X(t;6)dt. 

At Jt a +(n-\)M 



(4.14) 



The bin count x„ is a Poisson random variable whose probability is presented by 

p(x n = t,e)= ( n( d )Af ) exp(-A w (9)Af), fc = 0,l,2,... (4.15) 

and its mean and variance are 



E[x„;6] = var(x„;0) 
= Xn(8)At. 



(4.16) 



As jc„ random variables are independent, their joint probability mass function 
(pmf) equals 

N 

P {x-e)=Y[ P {xn,e) 



n=\ 



n«.,(-»). 



(4.17) 



«=i 



To derive the CRLB using (4.17), first, it must be checked if the regularity con- 
dition (4.5) holds. Therefore, the natural logarithm of (4.17) must be calculated 



lnp(x;9)= X (x„ln[A„(0)Af]-A„(0)Af-ln(x„!)). (4.18) 

H=0 

The derivative of (4.18) with respect to equals 

^M«*) = X (*-^'^W^-^(^) (4-19) 
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and its expectation is given by 

d 
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do 



lnp(x;0) 



N-l 



lX E[Xn] '^m'^ um ~^ l,{6]A ' 



de' 



N-l 



lX um -^-^ um -^ l,{6]A ' 



de' 



o. 



(4.20) 



Hence, the regularity conditions (4.5) hold. 

To obtain the CRLB, the Fisher matrix 1(9) must be found. From (4.7) and (4. 16), 
the matrix elements, /y, are given by 



I,j = -E 



N-l 



ddiddj 



lnp(x;0) 



,?,^ l " ll 's4* (8)i ' 1+ s4 y(,)i ') 



%{-^Hxm-me l K{m 



X n (6)At 



^U0)A t -^UB)A t } + J ^U0)A t ) 



N-l 



ZiJM-hiM-kM-*- 



(4.21) 



By letting At — ► 0, the summation in (4.21) can be replaced by an integral 



'/ de, 



3 X(t;8)~X(t;8) 



d6i 



X(t;6) 
Using (4.8), the partial derivatives in (4.22) are 



-6t, i,7 = l,2. 



(4.22) 



■^-X(f,e) = ^-X(f,e)=X s -h , (<j )Q + (t-t )f o ) (4.23) 

ddi d% 

^-X(t;6) = ^-X(t;9) = X s - (t -t ) ■h'(cj> + (t -t )f o ). (4.24) 

d0 2 df 
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Substituting (4.23) into (4.22), In may be found, 



/n = 



'f 



dipo 



A(r;0) 



-df 



/ 



Ht;6) 

'0+T obs [A s /z'((^o + (f-fo)/o)] 2 



A b + A s /!(0 o + (?~?o)/o) 



df 



T,** [X s h'(<h+fot)] 2 



-df. 



(4.25) 



The observation time can decomposed as r o b s =N„-P+T p , where P is the pulsar 
period, and < T p < P. Hence, the integral in (4.25) can be split as follows 



V » / mP [Wifo+fot)] 2 



£?l\J(n-i)Pfa+JLih(<pQ+f t) 



At 



NpP+Tp [Xsti(<$>0+f o t)\ 



N p P A b +Asft(0O+/ o f) 



At. 



(4.26) 



If the observation time is long enough, i.e., r o b s ^S> then r o b s ~ N p • P and the 
last integral in (4.26) almost equals to zero. Also, because h(-) is periodic, all of the 
remained integrals in (4.26) are equal. Therefore, 



p rnP 



hi 



[hti{fa + f t)} 2 

^1 J(n-l)P A b + X % h(<j)0+f o t) 



df 



■■Nr. 



[Kh'(<h+fot)f 



■At. 



p k %b + Xji{<fo+f t 
By changing the integration variable, and noting that f = \jP and r o b s « N„P 



(4.27) 



/ J<j> A b + A s /f(>) 



Tob 



Ab+Khtt) 



d(j> 



'obs 



[hti 



o A b + A s /i(0) 



A<j>. 



(4.28) 
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To find I>zi, (4.24) is substituted into (4.7), the integration variable is changed, 
and the integral is split into summation of the integrals in each period 



hi = 



'O + Tobs 



df 



Ht-,e) 



-At 



Ht-.e) 

'o+r obs [A s (r-f )/ i , ((/)o + (f-f )/ )] 2 
A b + A s /z((/)o + (r-fo)/o) 



At 



(4.29) 



T °b S [Is-t-tifa+fpt)] 

o A b + A s /z(0o + /of) 



-df 



n=\ 



nP 



L 



[h-t-h'{^+fot)\ 



'n,,p fa + L,h(<j)o + f t) 



At 



At. 



(4.30) 



As the observation time is long enough, the last integral in (4.30) is negligible. 
Also, notice that if (n — \)P < t < nP, then t in (4.30) can be replaced by 



t=(n-l + <j))P, 



(4.31) 



where < < 1 . Hence, 



/22 Sr [X s -t-h>(<j>o + f t)] 2 dt 
~~~ £J(."-i)p h + k s h(<t>o + f t) 



V P A 



[V(fI-l + 0)/»-Jl'(flD + (w-l + 0)/oP) 

A b + A s /!(0o + («-l + 0)/oP) 



-Pd0 



= **1 



'> yi(»-l + 0) 2 -[^A'(^, + »-l + 0)] : 



_,Jo 



B=l 



Ab + A s /;(0o + «- 1 



d</» 



= p3 | /•i (n-l + 0) 2 -[A g ^ + 0)] 2 
~i io A b + A s /? (0o + 0) 

By changing the integration variable in (4.32), 

3 ^ ^o+i(„-i + 0-0 o ) 2 .[A s /z'(0)] 2 



(4.32) 



«=1 



A b + A s A(0) 



df 



(4.33) 
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Because N p ^> 1, the integral value of the term that contains (n — 1 ) 2 is dominant. 
Therefore, 722 can be approximated by 

722 w p I / A I IdV d ^ (434) 

^ifc A b + A s /z((|>) 

which equals 

fa .^ W-l)W) W -l) /*« JWMla, (4 .35, 

and can be approximated by 

JV* *>+i [W0)] 2 

/22 ~ P T4 h+UW # 



4 S /■' [V*'W] 2 



3 Jo A b + As/i(0) 



(4.36) 



Substituting (4.23) and (4.24) into (4.7), changing the integration variable and 
splitting the integral, give rise to 

,.- ^-A(f;0)~A(f;0) 

/l2 = /21 = i We) dt 



I 



'o+r obs [W (fr + (t- fo)/o) ] [A s (r - f )/*' (ft + (t -to)f )] 



h + ^h(0Q + {t-t o )fo) 



o Ab + A s /z((^o + /oO 
|/^ f -[A^(0o + / o Q] 2 df 

„~1 \J(n-l)P h + k s h((j)o + fot) 

NpP+T P t-\l s h'(<bo + fot)] 2 



w p /> A b + A s /i((/)o + /oO 
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As N p ;§> 1, the last integral in (4.37) can be neglected compared to the first term. 
Using the change of variable (4.31), (4.37) can be simplified as 

«=i •/ (»- 1)^ Ab + ^h (<j>o + f t) 

* rl (n-l + $)P-[^h'(fo + (n-l + $)f P)] 2 
iJo h + hsh($ + (n-l + $)f o P) 

,;Ti Jo h + V*(0O + n - 1 + 0) 

^7o A b + A s /!(0 o + 0) 

°"ir" "T,g )f ^ <« 8 > 

„=l^o A b + A s /i(0) 

The term containing (« — 1) has the largest contribution in the summation (4.38), 
because A^>1. Therefore, (4.38) can be approximated by 

„=i^o A b + A s /!(0) 



=p' y ' ; i :y' d». (4.39) 

Since the term containing iV? is dominant, (4.39) can be approximated by 

in=hi~p 2 ^- ; 1 ; ;; J . s d< ^ (4 - 40) 

2 7^o A b + A s /z(0) 
and using the approximation r o b s « A/pP, it can be simplified to 

Z^^riMML^. (4 .41) 

2 7o Ab + A s «(^) 

From (4.28), (4.36), and (4.41), the Fisher matrix and its inverse can be found. 
These matrices are given in (4. 10) and (4. 12), respectively. ■ 
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If the observed frequency, /<,, is known then the only unknown parameter to be 
estimated is <po- In this case, from (4.28), the CRLB for estimation of (po is given by 

CRLB(fo) = — !— , (4.42) 

where I p is given in (4. 11). 

The pulse delay estimate, %, can be determined from (4.3) as 

?d = ^A (4.43) 
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Hence, if velocities of the spacecraft are known, fa does not need to be esti- 
mated, and from (4.42) and (5.1), the CRLB for estimation of t^ is given by the 
following [44]. 

2/1 



4.4 Discussion 



^-jlsi, 1 liM) 



Pulse delay estimator (4.43) is just a function of the velocity of spacecraft farther to 
the X-ray source. This may seem counter intuitive as one expects it to be a function 
of velocities of the both spacecraft. However, this is not the case. The pulsar signal 
is first detected on the closer vehicle. Then, assuming the clocks on both detectors 
are synchronized, since the signal is emitted by the X-ray source, the time duration 
it takes for the signal to be detected on the farther vehicle is just a function of its 
own velocity. 

Equation (4.44) shows that the CRLB value for estimation of the pulse delay is 
inversely proportional to the observation time. In other words, to obtain more accu- 
rate pulse delay estimates, longer observation time is needed. Furthermore, since the 
CRLB is inversely proportional to I p , which, in turn, is proportional to A s , making 
measurements on brighter pulsars results in more accurate pulse delay estimates. 
All of these analytical results are also consistent with intuition. 



4.5 Numerical Examples 

The CRLB for estimation of the relative distance is calculated for eight different 
X-ray pulsars. It is given by a — c^/CRLB(fd), where c is the speed of light. The 
employed pulsars are presented in Table 4.1, and their profiles in two cycles are 
shown in Fig. 4.2. The spacecraft velocities are assumed to be vj = 3 km/s and 
V2 = 9 km/s; and the initial phase observed at the first detector is (pi — 0.2 cycle. 
The relative distance between the spacecraft is assumed to be 180 km. Hence, the 
pulse delay is t^ — 0.6 ms. 



4.5 Numerical Examples 
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Table 4.1 Accuracy of the 
X-ray pulsar measurements 



Pulsar 



Period (s) 



min o~ (m) 



B0531+21 


0.0335 


2.6497E3 


B0540-69 


0.0504 


4.4901E3 


B0833-45 


0.0893 


1.8775E3 


B1509-58 


0.1502 


2.8866E3 


B1821-24 


0.0031 


188.68E0 


B1937+21 


0.0016 


27.5 16E0 


B1055-52 


0.1971 


1.1775E4 


J0437-47 


0.0057 


461.08E0 
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Fig. 4.2 Pulsars' profiles, h((j>), over two cycles 
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The observation time is assumed to be r o b s = 100 s, and a is calculated for 
different values of X\, and A s . The results are plotted in Fig. 4.3. As these figures 
show, the rate of change of the a surface drops for large values of A s . Note that 
from (4.44), the a value is inversely proportional to A s and the observation time, 
r o bs- Hence, to obtain smaller values of a for relative distance estimation, either the 
observation time or A s must be increased. The minimum values of the a surfaces 
are also given in Table 4. 1 . The minimum values change in a range of a few tens of 
meters to a few tens of kilometers. 
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Fig. 4.3 The a surface for r obs = 100 s 
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4.6 Summary 

In this chapter, we formulate the pulse delay estimation problem. We show that the 
pulse delay is composed of the true time delay between the detected pulsar rate 
functions (which is due to the geometric range) and the differential time between 
the detectors' asynchronized clocks. We present the CRLB for estimation of the 
spacecraft absolute velocities and the initial phase of the pulsar intensity functions. 
We also present the CRLB for estimation of the pulse delay when the spacecraft 
velocities are known. Some numerical examples are given as well. 



Chapter 5 

Pulse Delay Estimation Using Epoch Folding 



5.1 Introduction 

In this chapter we offer our first approach for estimation of the pulse delay. The 
proposed approach is to retrieve the photon intensity functions on each detector via 
epoch folding to obtain the empirical rate functions, i.e., Xk{t). Then, the empirical 
intensities will be used for estimation of the initial phase on each detector. Finally, 
from (4.3), the pulse delay can be estimated using the estimates of (j)\ and $2, 

? d = ^. (5.1) 

h 

Based on epoch folding, two different pulse delay estimators are presented in 
this chapter. Section 5.2 suggests using the correlation function between the em- 
pirical rate function and the true one to estimate the initial phase of the intensity 
function. In Sect. 5.3, a least-squares approach is proposed for estimation of the ini- 
tial phase. Some clarifying remarks are given in Sect. 5.4. The effect of absolute 
velocity errors on the performance of the proposed pulse delay estimators is studied 
in Sect. 5.5. We study the numerical implementation of the estimators in Sect. 5.6. 
Finally, we provide several numerical examples to assess the performance of the 
proposed estimators in Sect. 5.7. 



5.2 Cross Correlation Technique 

Let Tj, be the bin size, and Xj{ti) be the empirical rate function on y'th detector, where 
< tt < P. The known rate function is 

A(f,-) = A b + A s / ! ((f,--fo)/o). (5.2) 

Let the delayed rate function be 

X{t r ^j) = h + Kh{<$>j + {ti-h)f ). (5.3) 
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Therefore, 



where 
and 



X j (t i )=X{t i ;<S>j)+h j (t i ), .7 = 1,2, (5.4) 

E[nj{ti)}=0 (5.5) 



var[fy(f f )] = of (tt) 

= P L Wt,;4>j). (5.6) 

-*obs 

The initial phase <pj can be estimated by identification of the maximum of the 
cross correlation (CC) function between X(t) and Xdt) 

</>, = - arg max ^(i//), 7 = 1,2, (5.7) 

VG(0,1) 

where 

Rj{w)= T X*{t)lj{f,w)6t (5.8) 

is the CC function [44]. 

The asymptotic performance of this estimator is presented by the following 
theorem [44]. 

Theorem 5.1. The pulse delay estimator obtained from (5. 7) and (5.1), is unbiased 
for r o b s S> 0, and its asymptotic variance is given by, 

2 [ (X h + X s h(<j)))[X s h / ^)} 2 d<j) 

var[t A ] = -^ — ^— , (5.9) 

/ 2 2 r obs (J [Xsh'^fdij) 

where h'((j)) = dh/d(j). This means that the estimator is consistent since its variance 
converges to zero as more TOA measurements are incorporated. Furthermore, it 
is not asymptotically efficient, meaning that its variance does not converge to the 
CLRB infinite time, i.e., 

var[i d }> CRLB(t A ) for T ohs <c°. (5.10) 



Proof. First, the performance of the cross correlation initial phase estimator on each 
detector is analyzed. For simplicity in notations, the subindex /' of (j)j in (5.4) is 
dropped. Let 

Xi (t) = X{t) 

= X b +X s h((t-t )fa) (5.11) 
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and 

xi(t)=X(M). (5.12) 

Assuming <j> — —f D, x 2 {t) can be rewritten as 

x 2 (t) = A b + X s h((t - to) f + </>) + n(t) 
= h + AJi((t -t - D)f ) + h(t) 
= xi(t-D) + n(t). (5.13) 

Let D be the estimate of D. As the signals x\ (?) and x%{t) are band-limited, the 
cross correlation (5.8) has a derivative at D, and it follows from (5.7) that 

R(D)=0. (5.14) 

By linearizing R(D) about the true time delay D, the following is obtained [53]. 

R(D) + {D - D)R{D) = (5.15) 

Defining the delay estimation error 

S=D-D (5.16) 

from (5.16) and (5.15), it can be expressed as 

5 = -®. (5.17) 

R(D) 

The signals x\ (t) and x 2 (t) have the following spectra 

Si(f)=S x (f) (5.18) 

and, 

\s x [f)e-t 2 *f D + T b Y i 1i{k)<r jflnm l/l </s/2 
S 2 {f) — < ft (5.19) 

I |/| > / s /2, 

where A: = £#, / s = 1/7J,, and S^(f) is the spectrum of A(f), which equals 

Sxif) = 7iXA(^)e-^ 27r /^. (5.20) 

k 

The cross-correlation function is the inverse fourier transform of the cross power 
spectrum 

R(r)= f" Sl(f)S 2 (f)e j2nfr df. (5.21) 
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Therefore, 

*(t) = j2nfjS\{f)S 2 {f)J 2n f T df 

(5.22) 

Assuming h(k) is small (or equivalently T \, s ^> 0), (5.22) can be linearized about 
the true values of the samples as follows [53]. 

*M~I*(*)§§§. (5-23) 

where 

dR(x) _ .,,__„ ff*/ 2 
dn(k) 

Therefore, 



jlnnf* fSl(f)e-J 2 *f kT »eJ 2 Kf T df. (5.24) 

■/-A/2 



*(T)«./2*2i2*(Jt)( [~ fSiifloPtfWb-^A . (5 .25) 



Hence, 

= T h J j h{k)X{k). (5.26) 

From (5.21), 

*(t) - -4^_£] fS t l (f)S 2 (f)^ 2lc f' t df 

« -AK 2 f/\S x {f)\ 2 J 2 *f^- D )df. (5.27) 

Then, 

^(D)«-4^r/ 2 |5 A (/)| 2 df 



7*£A 2 (fe). (5.28) 
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Substituting (5.26) and (5.28) into (5.17), the estimation error is given by 

X«(*)A(*) 

8^^=^- ■ (5.29) 

EA 2 (£) 

k 

Therefore, from (5.29) 

E[8]=0. (5.30) 

Since <j> = —f D, this means that the phase estimator is unbiased as well, 

E[$]=j. (5.31) 

Furthermore, 



Xvar[«(fe)]-A 2 W 
2*W 



var[5] = ^— -=— , (5.32) 



where var[n(fc)] is given in (3.42). Therefore, 



var[<5] = =■. (5.33) 



Let T h — > (or equivalently Nb — ► °°), then (5.33) becomes 

1 rP 



var 



[8]: 



1 -?i> -/o 

Two 
j / A(f)A 2 (r)df 

' ' / A 2 (f)d? 



(5.34) 



The phase value varies in the range < <j> < 1, and the time span is < t < P. 
As a result, they are related as (j) = f/P. Therefore, 

dr = Pd0 (5.35) 
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and 

j t m = ^m- (5.36) 

Thus, (5.34) can be simplified as, 

I / A(0)A 2 (<^)d</> 
var[5] = ±- ■ ^° 

p y'wm 

p j\(<i>)k 2 ($)d$ 

Np (j[ I A 2 (0)d0 
f2 f A(0)A 2 (0)d0 

Jo JO 



(5.37) 



Since = —f Q D, 



l 

A(</>)A 2 O)d0 



/ o A 2 (<^)d(^ 

Therefore, substituting k(<j>) — Ab + A s /i(0) in (5.38), and using (5.1), the vari- 
ance of the time delay estimator is determined as given in (5.9). 

To compare the performance against the CRLB it is shown that var[fd]/CRLB(?d) 
is greater than one. Using (5.9) and (4.44), the ratio equals 



var[f d ] '- SiW^- L g3(0)d(/> 



CRLB(f d ) 



/ g20)# 

Jo 



(5.39) 



where 



gi(0) 4 (A b + A s / J (0))[A s / ! '(</.)] 2 (5.40) 

g 2 (tf>) 4 [A s /*'(</>)] 2 (5.41) 

»w = AbTi^y (5 - 42) 
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Note that from Cauchy-Schwartz inequality, 



l^g^d^-^g^d^y^^/g^g^d^j 



(5.43) 



From (5.40), (5.41), and (5.42), it can easily be seen that 

gift) = V8i(<t>)g3(<t>)- (5-44) 

Therefore, (5.43) and (5.39) result in 

var[? d ]>CRLB(f d ). (5.45) 

The equality holds if and only ii g^((j>) is a multiple of g\ (0), which implies 

/Lb + A s /z(0) = const. (5.46) 

Since this is not the case for a general h(<p), equality does not hold, and the 
variance does not attain the CRLB asymptotically. ■ 

5.3 Nonlinear Least Squares Technique 

The initial phase can be estimated by fitting the empirical rate function, A/(f), to 
the true known rate function, k(t), via solving a nonlinear least squares (NLS) op- 
timization problem [43]. The NLS cost function is defined based on minimizing the 
difference between the empirical rate function and the true one at the j'th detector, 
i.e., Xj(t) — X(t;<j)j). Hence, the objective function, J(<j>j), is defined as follows. 

J(<l>j) = 1 LMi)-Mti;<l>i)) , 7 = 1,2 (5.47) 

The unknown initial phase is estimated by minimizing the cost function (5.47), as 

7 = argmin/(0 ; ), j=l,2- (5.48) 

^e(o.i) 

The following theorem summarizes the performance of the proposed pulse delay 
estimator [44]. 

Theorem 5.2. The pulse delay estimator obtained from (5.48) and (5.1) is asymp- 
totically unbiased, and its asymptotic variance is given by (5.9). That is, the NLS 
estimator's asymptotic performance is the same as the CC estimator. 
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Proof. For simplicity in notations, the subindex j in (5.47) is dropped and the cost 
function is restated as, 

/(0) = (x-f(0)) T (x-f(0)), (5.49) 

where 



and 



x^{X( tl )l(t 2 )---l(t Nb )) T (5.50) 



f(0)^(A(r i; 0)A(f 2 ;</.)---A(f Wb ;</.)) T . (5.51) 



Now, assuming that the selected minimum is in the correct neighborhood (or 
equivalently the observation time is long enough), the mean and variance of the 
estimator can be approximated by linearizing f(0) about the true value of phase, 
0°, as 

f(0)«f(0°) + (0-0°)h, (5.52) 



where 

d(j) 

Therefore, (5.49) becomes 



(5.53) 



7(0) = (x - f(0°) + h0° - h0) (x - f(0°) + h0° - h0) . (5.54) 

Noting that x — f(0°) + h0° is fixed and is the variable, the linear least-squares 
solution to the cost function (5.54) is given by 

0=yt- 1 h T (x-f(0°) + h0°) 
= 0° + r- 1 h T (x-f(0°)), (5.55) 



where 



r = h T h. (5.56) 



The estimation error is defined as 

^ = 0-0°. (5.57) 

Hence, from (5.55) 

E[e+] = r- 1 h T (0°)£[x-f(0 )] (5.58) 

and 

var[e ] = r _1 h T var[x-f(0°)]hr _1 . (5.59) 

Note that from (3.40) 

x-f(0°)=n, (5.60) 

where 

n= {n{h)h(t 2 ) ••• «(fiv„)) • (5.61) 
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Hence, employing (3.41) and (3.42), the asymptotic mean and variance of the 
estimation error, as r o b s — ► °°, are 



and 



where 



£[«*]= (5.62) 



var[e^} = P^r- l h T Lhr-\ (5.63) 

*obs 



L4diag(A(r i; (/. o ),A(r 2 ;(/) o ),---^(^;0 o ))- (5-64) 



This means that the NLS estimator % is asymptotically unbiased. By substitut- 
ing the values of r, h, and L into (5.63), the asymptotic variance of initial phase 
estimator is given by 



A'l 



var[&] = — ■— =-. (5.65) 

(£(A'(r, ;( n) 2 J 

By letting N\, — ► °° and substituting the value of A(-) from (3.27), (5.65) can be 
replaced by 

/ (A b + A s /z(</.))[A s / i '(0)] 2 d<^ 
\ar[$j] = i5 _^_. (5.66) 

r obs H [AsA'Wfd^. 

Therefore, from (5.1), the time delay estimator's asymptotic variance equals 
2var[0j]//2 , which is given in (5.9). ■ 



5.4 Discussion 

In [54], it is stated that if some regularity conditions are satisfied, the NLS estimator 
is consistent. It is easy to check that the mentioned conditions are met by the NLS es- 
timator (5.48). Furthermore, it is mentioned that the NLS estimator is not in general 
asymptotically efficient, and this property depends on the underlying distribution of 
the noise. The significance of Theorem 5.2 is that it states (5.48) is indeed not an 
asymptotically efficient estimator. A meaningful measure to compare the estimator's 
performance to the CRLB is the asymptotic relative efficiency (ARE), which is the 
ratio between the estimator's asymptotic variance and the CRLB [55]. The ARE for 
the proposed NLS estimator is calculated and given in (5.39). From (5.39), it is clear 
that the LNS estimator's ARE is only a function of the pulsar profile, h(<j>), and the 
photon constant rates, Ab and A s . 
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Also note that the proposed NLS phase estimation approach is different from the 
one given in [56], which is typically used for radio pulsar timing. The method in 
[56] is based on minimizing the difference between the Fourier transforms of the 
empirical rate function and the true one. 



5.5 Absolute Velocity Errors 

The presented pulse delay estimators were analyzed based on accurate knowledge of 
the spacecraft absolute velocities. If there are absolute velocity errors, the proposed 
pulse delay estimation approach still works under certain conditions. If the velocity 
errors in the direction of the X-ray source are small enough, the resulted pulse delay 
errors can be neglected. This is a result of the upper bound presented in (3.72). 
If the velocity errors on each spacecraft are almost equal, Avi ss Av2, then from 
(3.66), A0i rs A02- This results in the same deterioration of initial phase estimates 
on each detector. Hence, from (5.1), as the difference between the phase estimates 
is employed for estimation of the pulse delay, it remains asymptotically unbiased. 
However, the variance of estimation error becomes larger than the case where the 
absolute velocities are perfectly known due to two reasons [44]. One is the effect of 
the velocity error on fa in the denominator of (5.1). It is given by 

A<S 
Af d = --fA/ 2 , (5.67) 

h 
where 

A0 = 0i - <fc (5.68) 

and 

Av? 
A/ 2 =/s— . (5.69) 

c 

Therefore, the time delay estimate error due to this velocity error equals 

"AV2 

Af d = -A0— =•/„ 
c fi 

» Avo 
«-A0 — ±. (5.70) 

cfa 

This causes cAfj meters of error in position estimation 

8 Ax = cAtd 

--A0-A (5.71) 
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The second reason that the variance of the pulse delay estimator increases is that 
the variance of the initial phase estimates becomes larger. These estimation errors 
can be modeled by increasing the variance of the measurement noise in the Kalman 
filtering stage. 

If the aforementioned assumptions about spacecraft velocities are not valid, the 
absolute velocities must be estimated as well as the initial phase values. We propose 
two solutions for this problem. 

One approach is to simultaneously estimate /# and fa. This is possible by 
using the TOA pdfs (6.1) to form two likelihood functions, and then solve the 
two-dimensional ML estimation problems. More details about this approach are dis- 
cussed in Chap. 6. 

In a different approach, we propose to employ the geometry of pulsars in the 
sky to estimate the spacecraft absolute velocities. For more details on this method, 
please see Chap. 7. 



5.6 Computational Complexity Analysis 

5.6.1 Epoch Folding 

To use the CC or NLS initial phase estimators, the empirical rate function must be 
obtained through epoch folding. To obtain it, from (3.39), N\,(N P — 1) additions, and 
2N\, divisions are needed. Hence, the total number of calculations for epoch folding 
is N\, (N p + 1 ) flops. This shows that the amount of necessary calculations for epoch 
folding is a linear function of N p . In other words, it roughly linearly grows as the 
observation time becomes longer. 

5. 6. 2 CC Estimator 

As X(t;y/) is a discrete signal which is available at multiples of Tf, s, the discrete 
cross-correlation function must be calculated. It is given by 

1 N b 
MW) = w Y. x ^ kT b)x2(kT b ;\ i r) 1 (5.72) 

/Vb k=l 

where X\ (?) and xjit) are defined in (5.11) and (5.12), respectively [57]. 

The cross-correlation function can be computed in the time domain using (5.72), 
or in the Fourier domain 

Rd(w) = &- l {X\{fk)Xl{fk)}> (5-73) 

where X\ (//.) and ^(/j) are the discrete Fourier transforms of x\ [kTj,) and JC2(Wi). 
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If the resolution of the initial phase estimate is needed to be finer than 1/iVb, 
it is necessary to interpolate the cross-correlation function. One approach is to ap- 
proximate the correlation function by a convex parabola in the neighborhood of its 
maximum, as 

R D (\lf)=a\i/ 2 + by + c, (5.74) 

where a, b, and c are parameters fitting the measured correlation. Using this approx- 
imation, the continuous time-delay estimation can be performed by finding the apex 
of the parabola, 

b, = ~. (5.75) 

la 

Hence, to obtain the continuous initial phase estimate, two steps are needed: 

• Locate the index k m of the maximizer delay Ro(k m T s ), where T s = 1/N\,, to find 
the coarse estimate. 

• To find the subsample estimate, use the maximum cross-correlation delay 
Ro{k m T s ), and two of its adjacent values Ro(k m T s — T s ) and Ro(k m T s + T s ) 
for a parabolic approximation, and add it to the coarse estimate, 

* = . T _ 1 R D {k m T, + T,)-R D {k m T,-T,) 

9 ms 2' R D (k m T s + T s )-2R D (k m T s )+R D (k m T s -T s y 

Note that due to the parabolic approximation, the phase estimate (5.76) becomes 
biased [57]. 

Using (5.72), to obtain the discrete cross-correlation function, for a fixed I//, 
N\, multiplications, N\, — 1 additions, and 1 division are performed. Since the 
search domain to maximize Rd(w) i s (0; 1) cycle* which is divided into iVg grids, 
a total number of 2N g N\, flops must be performed to obtain the initial phase 
estimate. 

From (5.76), we can see that three additions and two multiplications are needed 
for interpolation. It is also clear that the interpolation is performed only to find the 
maximizer point. Hence, the number of calculation flops is not a function of the 
observation time, r o b s - 

In summary, to obtain the CC estimator, first epoch folding is performed which 
needs N\,(N p + 1) operation flops. Then the search interval (0, 1) cycle is divided 
into Ng grids, and the correlation function is obtained for each grid. This imposes 
2N\,Ng operations. 

If the correlation function is obtained using the Fourier transform from (5.73), 
first the Fourier transform of each signal must be found. Then N\, multiplications are 
performed and using inverse Fourier transform on the resulted signal, the correlation 
function is determined. As fast Fourier transform techniques are widely available, 
this approach is expected to be faster than the time domain method. 
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5. 6. 3 NLS Estimator 

To construct the NLS cost function, using (5.47), N\, subtractions, N\, multiplica- 
tions, and N\, additions must be performed for a fixed 0o. Hence, the total number of 
floating-point operations is 3A^- This means that the computational cost to construct 
the NLS cost function is not a function of the observation time. 

Although for any given 0o, A (f,, (j>o) is known, the mathematical equation of the 
rate function A (t ; 0o) is not usually available in practice. Hence, the gradient of the 
cost function, if needed, must be calculated numerically. Furthermore, the cost func- 
tion is not convex in general. This means that it usually has multiple minima. Hence, 
to avoid getting trapped in local extrema, and since (5.48) is a scalar optimization 
problem, employing a one-dimensional grid search over the interval (0, 1) cycle is 
necessary to solve it. In other words, first, the epoch folding procedure, with a cost 
ofNb(N p + 1) flops, is performed. Then, the search interval (0, 1) is divided into N g 
grids. For each grid, the cost function (5.47) is calculated, and the minimizer grid is 
found. This imposes a total floating-point operations of 3NgNb flops. 



5.7 Numerical Examples 

First, performance of the initial phase estimators on one detector is examined. 
Then, a relative navigation scenario is presented and the pulse delay estimators are 
simulated. 



5. 7. 1 Initial Phase Estimators 

The Crab pulsar (PSR B0531+21) is chosen as the X-ray source. The constant ar- 
rival rates are chosen to be Ab = 5 and A s = 15 ph/s. The spacecraft velocity on the 
direction vector pointing to the pulsar is assumed to be v = 3 km/s; and the initial 
phase observed at the detector is (pi = 0.2 cycle. A Monte Carlo simulation was per- 
formed in which 500 independent realization of photon TOAs were processed by 
the phase estimators. 

To evaluate each estimator's performance, the empirical rate function is derived 
for each observation time, and the initial phase is estimated by optimizing the cost 
functions given in (5.8) and (5.47). The optimization is done using a grid search 
algorithm over [0,1) cycle. The interval is divided into 1,024 grids and the cost 
functions are calculated for each one. Note that the magnitude of the estimation 
error is calculated modulo one cycle, i.e., 

|e|=min{mod(0o-$o,l), mod (<fo- 0o,l)}- (5.77) 

For instance, if 0o = 0. 1 and </>o = 0.9, the error is 0.2 cycle, and not 0.8 cycle. 
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Fig. 5.1 RMS error for the CC and NLS phase estimators 



In Fig. 5.1, the root mean square (RMS) of the CC and NLS estimation errors for 
each observation time are plotted against the square root of the CRLB(0o) and the 
analytical value, presented in (4.42) and (5.66), respectively. The mean square error 
(MSE) of the estimator which measures the average mean squared deviation of the 
estimator from the true value, is defined as follows. 



MSE(0o)=£[(<£o-4>o) 2 ] 



(5.78) 



The MSE can be rewritten as 



MSE(<fc>)=£<^ U)-£($ 

= var(0 o ) + [E(^q)-(J)q\' 
= var(<£ o ) + /? 2 (0o), 



E($o) 



(5.79) 



where 



b{<k)) =#(&)- 0o 



(5.80) 



is the bias of the estimator. This shows that the MSE is composed of errors due to 
the variance of the estimator as well as the bias. 

The RMS plots in Fig. 5.1 show that as time goes on, the variance of estimation 
error approaches zero, implying that both estimators are consistent. The difference 
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between the RMS error and the CRLB is determined by the difference between 
(5.66) and (4.42). It is in general a function of the pulsar rate function, but it is 
always inversely proportional to the observation time. As the observation time is 
reduced, there is a threshold point at which the RMS error starts to deviate from the 
CRLB value. The reason is that as the observation time goes below this threshold, 
the realization of TOAs does not represent the rate function precisely. Therefore, 
they result in a distorted cost function whose extremum happens at a farther point 
to the true parameter, (j)". As a result, the estimator becomes biased, and as (5.79) 
shows, the estimation error's variance starts to depart from the CRLB. 

The CC cost function (5.8), obtained by Monte Carlo simulation over 500 re- 
alizations of TOAs, is plotted in Fig. 5.2 for different observation times. The plots 
show that for short observation times, the CC function has a smaller peak, and as the 
observation time increases and more TOAs are measured, the peak value becomes 
larger. 

The cost function (5.47) is also plotted in Fig. 5.3 for different observation times. 
The plots show that as the observation time increases, the cost function becomes 
smaller. 

Figures 5.2 and 5.3 also show that each of the cost functions has multiple 
extrema. Hence, they verify that, to avoid getting trapped in local extrema, adopting 
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Fig. 5.2 The CC cost function 



64 




5 Pulse Delay Estimation Using Epoch Folding 
x10 6 



3.155 




x10 




x10 



0.5 

0o (cycle) 

Fig. 5.3 The NLS cost function 
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a grid search optimization approach is necessary. As the search is performed over 
the whole interval (0, 1) cycle, the algorithm initialization issues do not arise. 



5.7.2 Pulse Delay Estimators 



Now we simulate the pulse delay estimators. The spacecraft velocities on the direc- 
tion vector pointing to the pulsar are assumed to be vi = 3 km/s and vj = 9 km/s; 
and the initial phase observed at the first detector is <j)\ — 0.2 cycle. The Doppler 
frequencies depend on vi and V2, which are scalars. Hence, we do not need the 
spacecraft velocity vectors for this simulation. The relative distance between the 
spacecraft is assumed to be Ad — 180 km. Hence, the pulse delay to be estimated is 
fd = 0.6 ms. Since t<i is less than the pulsar period, there is no phase ambiguity. 

An observation time of r o b s = 100 s is selected to verify the epoch folding re- 
sults. For each observation time, photon time tags are folded back into one single 
pulsar cycle which is divided into N\, = 1 ,024 bins, and the rate functions on each 
detector are derived using (3.39). The empirical rate functions along with the true 
ones are plotted in Figs. 5.4 and 5.5. As the plots show, the empirical rate functions 
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Fig. 5.4 A 1 (f)forr obs = 100s 
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Fig. 5.5 X 2 {t) for r obs = 100 s 
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Fig. 5.6 The variance of A\ (t) for r o t, s = 100 s 



converge to the true functions. In Figs. 5.6 and 5 .7, the variance of n(t), obtained via 
simulation, is plotted and compared to the analytical variance for each detector. 

The effect of imprecise absolute velocity data on epoch folding is studied as well. 
The errors are intentionally chosen to be big: Avi = 1 ,000 m/s and Ai>2 = 1 ,250 m/s. 
For the selected parameters, the upper bound (3 .72) roughly equals 10 m/s. Since the 
velocity errors are bigger than the bound, the empirical rate functions are expected to 
be deteriorated. Figures 5.8 and 5.9 verify this fact. Comparing them to the true rate 
functions shows that they are shifted, and their peaks are flattened. The same phe- 
nomenon happens for the noise variance which can be seen in Figs. 5.10 and 5.11. 

Figures 5.12 and 5.13 show the RMS error of the proposed pulse delay esti- 
mators, obtained through simulation, against the CRLB and the analytical value, 
calculated as in (4.44) and (5.9). The mean square error (MSE) of the estimator is 
defined as 



MSE(? d ) =E[(k-t d y 



(5.81) 



and the cost functions are optimized using a grid search approach in the domain of 
(0,1) cycle. 

As expected, the RMS plots show that as time goes on, the variance of each pulse 
delay estimator approaches zero. Also, as the observation time is reduced, there is 
a threshold point at which the estimation error variance starts to deviate from the 
CRLB. It is due to the fact that the phase estimators become biased when there are 
not enough photon TOAs to process. 
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Fig. 5.7 The variance of h%(f) for r o t, s = 100 s 
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Fig. 5.8 Ai (0 when Av 3 = 1 ,000 m/s 



68 



5 Pulse Delay Estimation Using Epoch Folding 



150 



100 




15 20 

Time, t (msec) 



30 33 



Fig. 5.9 X 2 (t) when Av 2 = 1,250 m/s 
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Fig. 5.10 var[»(r)] when Avi = 1,000 m/s 
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Fig. 5.11 var[n(r)] when Av 2 = 1,250 m/s 
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Fig. 5.12 The RMS error of the NLS-based estimator 
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Fig. 5.13 The RMS error of the CC-based estimator 



Note that we performed this simulation to verify that the analytical RMS values 
match the numerical results. Different choices of A D , A s , vi, V2, <j>i, and Ad, only 
change numerical values of CRLB, analytical RMS errors, and simulation RMS 
errors. 

Finally, the performance of the pulse delay estimators is studies when the space- 
craft velocity errors are Avi = 1 ,000 m/s and Av2 = 1 ,250 m/s. The RMS error plots 
are presented in Figs. 5.14 and 5. 15. As the plots show, the asymptotic performance 
of the estimators degrades compared to the case where the absolute velocities are 
perfectly known. 



5.8 Summary 



We explain how to use the epoch folding procedure for estimation of the pulse delay. 
We presented two different initial phase estimators which both work based on epoch 
folding. The key idea in employing epoch folding for X-ray pulsar-based navigation 
is to recover the photon intensity functions on each detector and estimate their initial 
phase. One estimator works based on maximizing the cross correlation between the 
empirical rate function and the true one. The other one is obtained by minimizing the 
difference between the rate functions via solving a least-squares problem. We ana- 
lyze theses estimators and show that they are not asymptotically efficient but they 
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Fig. 5.14 The RMS error of the NLS-based estimator in presence of absolute velocity errors 
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Fig. 5.15 The RMS error of the CC-based estimator in presence of absolute velocity errors 
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are consistent. We investigate the numerical implementation aspects by analyzing 
the number of floating point operations. Further numerical examples on computa- 
tional complexity of the estimators will be presented in Chap. 6. We explain that 
the proposed pulse delay estimation technique works not only when the spacecraft 
velocities are perfectly known but also under some conditions where the velocity 
data is not perfect. By offering different simulation scenarios, we also examine the 
analytical results via numerical examples. 



Chapter 6 

Pulse Delay Estimation via Direct Use of TOAs 



6.1 Introduction 

An important point regarding the pulse delay estimation techniques introduced in 
Chap. 5 is that they need to employ the epoch folding procedure which needs the 
exact knowledge of the spacecraft velocities. Furthermore, it was shown that the pro- 
posed epoch folding-based estimators are not asymptotically efficient. To address 
these problems, we propose another method in this chapter. This approach is based 
on direct use of the measured TOA sets, {t\ },2\ and {t\ }/_?i • Because the mea- 
sured TOAs are being used directly, it is not necessary for the estimators to have 
access to the velocity data. Nonetheless, the effect of imprecise velocity data on 
the performance of the pulse delay estimator is studied. Computational complexity 
study is also performed. 

We formulate an ML problem in Sect. 6.2 for estimation of the Doppler frequen- 
cies and the initial phase values. Some remarks on numerical implementation of 
the proposed estimator are offered in Sect. 6.3. The ML computational complexity 
analysis is presented in Sect. 6.4. All computational complexity analysis results are 
summarized in Sect. 6.5. The effect of absolute velocity errors on the pulse delay 
estimator's performance is studied in Sect. 6.6. Numerical simulations are provided 
in Sect. 6.7. 



6.2 Maximum-Likelihood Estimator 

Employing the pdfs associated with the detected time tags, an ML estimation prob- 
lem can be formulated to estimate fa and fa . According to (3 . 1 0), the pdfs associated 
with each set of time tags are given by, 
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where 



Mk>fk)= [ f h{f,hJk)dt, 

Jtn 



£=1,2. (6.2) 



Recognizing the pdfs given in (6.1) as likelihood functions, the maximum- 
likelihood estimators (MLEs) are provided by maximizing each likelihood function 
with respect to the unknown parameters. Equivalently, the natural logarithm of the 
likelihood function or the log-likelihood function (LLF) can be maximized. 

M k 

LLF(fa,f k ) = J j ln(X k (t i ;<l) k ,f k ))-A(fa,f k ), *=1,2. (6.3) 

;=i 

If the observation time is long enough compared with the pulsar period then 
A[fa,f k ) in (6.3) shows a minimal dependence on the parameters fa and f k , that is, 
dA/dfa — 0, and dA/dfk — 0. The reason is as follows. Let the observation time 
contains N p pulsar cycles, T \, s = N p P+T p , where < T„ < 1. Using periodicity of 
X(t;fa,fk), (6.2) can be res-stated as follows. 



A(k,/*)= [° ' ' X{t;fa,f k )dt+N p [ X{t;fa,f k )dt. (6.4) 

Jtn JO 



Since X (■) > is a periodic function, its integral over one period is not a function 
of the initial phase and the Doppler frequency. As a result, the second term on the 
right hand side of (6.4) is independent of fa and fa, while the first term is a function 
of fa and fk- Hence, if N p 3> 1, then the second term is strongly dominant, and the 
value of A(0fc,//t) will not be sensitive to fa and/jt. In other words dA(fa)/dfa «0 
and dA (fk)/dfk ~ 0. Therefore, it can be dropped from the objective function, and 
the likelihood functions, denoted by *F(fa,fk), become, 

M k 

V{<$>k,fk) =J,ln(X k (t i ;faJ k )). (6.5) 

i=i 

The unknown parameters fa and f k can be determined by solving the following 
optimization problems, 

(fajk) =argmaxf(^,/j) ! k=l,2 (6.6) 

faJk 

and the pulse delay may be estimated from (4.43). 

The statistical properties of the MLEs (6.6) are summarized in the following 
theorem [51]. 

Theorem 6.1. The MLE determined from (6. 6) is unbiased for r o b s S> 0, and attains 
the CRLB, presented in (4.12). In other words, it is asymptotically efficient. 

Proof. To investigate the properties of the MLE, first note the following theorem 
[50]. 



6.2 Maximum-Likelihood Estimator 
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Theorem 6.2. If the pdf p(x;9) of the data x satisfies the regularity conditions 
given in (4.5) then the MLE of the parameter 6 is asymptotically unbiased and 
asymptotically attains the CRLB. It is therefore asymptotically efficient. 

Let 6 be the unknown parameter to be estimated. To apply Theorem 6.2, it must 
be shown that the regularity condition (4.5) is satisfied, where x = {?;}^i- From 
(3.10), 



/:/ 



lnp{x;6) = -A(9) + ^lnX{ti-,e). 



(6.7) 



Therefore, 



^M^) = 4A (e)+ lJ^.^(,9). (6.8) 



Because the observation time is assumed to be long enough, the first term in (6.8) 
vanishes and, 



do 



lnp(x;9) 



M 



i=\ 



1 



A(r,-;0) 36 



A(/,-;0) 



(6.9) 



To calculate the stochastic expectation appearing in the right-hand side of (6.9), 
the following corollary may be used [41]. 

Corollary 6.3. Let {tj}j =i be the photon TOAs with the pdf given in (3.10). Also, let 



.7=1 



(6.10) 



where r(t) is any general function. Then, the stochastic expectation of g ( {tj}j =i ) 
for a fixed N equals, 



,-A 



N\ 



E N [g({tj}U)]-—^ 
where A is defined in (3.11), and, 



(6.11) 



T^ / r(t)k{t)dt. 
Ji 



(6.12) 



Proof ( Proof of Corollary 6.3). Let Qpj be the event of receiving /V photons at any 
/V different increasing time instants in an interval [<o,f/]. Then, using the TOA pdf, 
given in (3.10), 



E N[g({tj}U)] = [ p({ti}li,N)flr(tj)dh-dt N , (6.13) 
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where Qn is the event of receiving N photons at any N different increasing time in- 
stants in interval [to , tA . Note that the sequence {ti , fe, •••;%} is in increasing order. 
As there exists A! permutations of ti, for a fixed A, the probability that A number of 
ti's occurs in no special order is N\ times that of the sequence occurring in increasing 
order. Therefore, 



A' 



p({ti}l u N)Y[r(t J )dt l ---dt h 



1 /•'/ /•'/ r'f 

Nil Jt '"' 

*—A rtf rtf rtf JV 



' p{{ti}UN)Y[r{tj)Atx 

'0 /=! 



Y1 , , , Y[r{tj)X{tj)At v 

A! Jt Q Jt J tg j= i 



■dt N 



^1 

A! 



r(t)X(t)6t 



r* 



■dt N 



(6.14) 

D 



Now it can be seen that the stochastic expectation in the right-hand side of (6.9) 
is a special case of (6.1 1) where N — 1. Hence, 



1 



k(tr,0) d9 



Kn-d) 



- A r 



r(t;6)X(t;8)dt, 



(6.15) 



where 

which yields in, 
E 



1 



A(r,-;0) d6 



r(t;9) 



Ku-,d) 



X(t;9) 

MtiWj 



,-A 



(i(t r ,e)-x(t ;e)). 



(6.16) 



(6.17) 



As A increases as a function of observation time, the term e~ A approaches zero, 
while X(tf;G) — X(to;0) remains finitely bounded in magnitude. Therefore, the 
right-hand side of (6.17) approaches zero. This means that the regularity condition 
(4.5) holds, and as a result, according to Theorem 6.2, the MLE is asymptotically 
efficient. ■ 
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Assuming the velocities of the vehicles are known, the Doppler frequencies do 
not need to be estimated, and the phase estimates may be obtained by solving the 
following optimization problem, 

$t = argmax *?(&)> k =1,2. (6.18) 

<fee(0,l) 

Furthermore, the pulse delay estimator is determined from (5.1). 

The statistical properties of the resulting pulse delay estimator are presented by 
the following theorem [49]. 

Theorem 6.4. The pulse delay estimator determined from (6.18) and (5.1) is 
asymptotically efficient, and its variance is given by (4.44). 

Proof. From Theorem 6.1, the phase estimators are asymptotically efficient. Hence, 
using (5.1), the pulse delay estimator is asymptotically efficient as well, and its 
variance equals 2CRLB($o) / '/| where CRLB(0o) is presented in (4.42). ■ 



6.3 Numerical Determination of the MLE 

To find the MLE of a parameter 6, the likelihood function needs to be maximized. 
If the allowable values of 6 lie in a certain interval then the maximum of likeli- 
hood function in that interval needs to be found. The safest way to do this is to 
perform a grid search over that interval. As long as the spacing between 6 values 
is small enough, it is guaranteed to reach to the maximum of likelihood function. 
However, if the range of 9 is not known or confined to a finite interval then a grid 
search may not be computationally feasible. In such situations, iterative optimiza- 
tion procedures are employed. A typical one is the Newton-Raphson method. In 
general, these methods result in the MLE if the initial guess is close to the true 
value. If not, they may not converge or only convergence to a local maximum is 
attained. An important point about MLE is that the likelihood function changes for 
each data set. 

The iterative methods attempt to maximize the likelihood function by finding the 
zero of its derivative, 

^ln/»(jr,0)=O. (6.19) 

The Newton-Raphson iteration is [50], 

(6.20) 



e^ = eV-(^ P (,,e)y^n P{x ;e: 



)=<#) 
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where 



d 



2 



dede T 



lnp(x;6) 



d 



2 



dOide 



\ap{x,e) (6.21) 



is the Hessian of the likelihood function, and d In p(x;9)/d6 is the gradient vector. 
In implementing (6.20), inversion of the Hessian matrix is not required because it 
can be rewritten as, 



d 



2 



30dO T 

d 2 



In p(x;0 

lnp(x;9 



g(*+i) 



dede 



e=ew d9 Fy ' 



( 6 - 22 > 



Regarding Newton-Raphson iterative procedure, one should notice that it may 
not converge. In particular, this may be the case if the Hessian matrix is small. Even 
if it converges, the point found may only be a local optimum. Hence, to avoid these 
possibilities, it is beneficial to use several starting points. 



6.4 ML Computational Complexity Analysis 

From (6.5), to construct the ML cost function for a fixed 0q, M additions must 
be performed. Also, the cost function must be evaluated at the measured TOAs. 
If the rate function's analytical equation is known, this does not impose a significant 
amount of computational cost. 

As the mathematical equation of the rate function is not usually available in 
practice, its value at the TOAs must be found by interpolation. If a simple linear 
interpolation is used, it costs four additions and two multiplications for each TOA. 
Hence, the total computational cost, for a fixed </>o, is 6M flops. Similar to the NLS 
or CC cases, because the cost function is not concave in general, adopting a grid 
search optimization approach over (0, 1) cycle is desirable. If the search interval is 
divided into N g grids, the total imposed computational cost due to the interpolation 
is 6MN g flops. This means that the amount of interpolation computations also grows 
as the number of received photons, M, increases. 



6.5 Computational Complexity: Summary 

Note that M ^> N p for long observations. Hence, compared to the NLS and CC 
cases, the amount of calculations for the ML estimator significantly increases as 
the observation time becomes longer. This is a noticeable disadvantage of the ML 
approach compared to the NLS estimation. This is to be avoided for implementation 
purposes. Hence, if the NLS or CC estimator's variances are within an acceptable 
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Table 6.1 Computational cost 








Estimator Additions 


Subtractions 


Multiplications 


Divisions 


CC N b {N p -l) + (N b -l)N g 
NLS N b (N p -l)+N h N g 
ML M + 4MN g 


N/A 

N h N g 

N/A 


N b N g 
N b N g 
2MNo 


2N h + l 

2N h 

N/A 
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distance from the CRLB, it is more convenient to implement one of them. This is 
indeed the case. It is clear from (4.44) and (5.9) that the difference between the NLS 
estimator's variance and the CRLB becomes smaller as the observation increases. 
As a result, especially for long observation times, implementation of epoch-folding- 
based algorithms is computationally more efficient. Note that the ARE, given in 
(5.39), can also be calculated for each pulsar to measure how far is the estimator's 
performance from the CRLB. 

If the cross-correlation function is calculated using the time domain approach, the 
computational cost for CC and NLS estimators is almost the same. As the correlation 
function can be found using the Fourier domain approach, its computational cost is 
expected to be less than the NLS estimator. 

All computational cost results, presented in Sects. 5.6.2, 5.6.3, and 6.4, are sum- 
marized in Table 6. 1 . 



6.6 Absolute Velocity Errors 

Similar to the CC and NLS cases, if the velocity errors on each spacecraft are almost 
equal, Avi « Av2, this results in the same deterioration of initial phase estimation on 
each detector. Hence, from (5.1), since the difference between the phase estimates 
is employed for estimation of the pulse delay, it remains asymptotically unbiased. 
However, the variance of estimation error becomes larger than the case where the 
absolute velocities are perfectly known. This error can be modeled by increasing the 
variance of the measurement noise in the Kalman filtering stage. If these assump- 
tions are not valid, the absolute velocities can be estimated from (6.6). 



6.7 Numerical Examples 

The algorithm presented in Sect. 3.8 is used to generate the photon TOAs associ- 
ated with the Crab pulsar (PSR B0531+21). The simulations are performed using 
the Monte-Carlo technique with over 500 independent realization of photon TOAs 
for each observation time. The constant arrival rates are chosen to be Ab = 5 and 
A s = 15 ph/s. 

First, the ML phase estimator is numerically simulated on one detector. Then the 
pulse delay performance is numerically examined. Finally, the computational cost 
needed for calculation of all of the proposed estimators is studied. 
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6 Pulse Delay Estimation via Direct Use of TOAs 



6. 7. 1 ML Phase Estimator 

The spacecraft velocity on the direction vector pointing to the pulsar is assumed to 
be v = 3 km/s; and the initial phase observed at the detector is 0i = 0.2 cycle. 

Figure 6.1 shows the RMS error of the MLE. The ML estimates are determined 
by solving the optimization problem (6.18). A grid search approach on the interval 
[0,1) cycle over 1,024 grids is utilized to solve (6.18). The ML plot shows that 
the estimator attains the CRLB for long observation times. Similar to the CC and 
NLS estimators, as the observation time drops below a certain threshold, the RMS 
error deviates from the CRLB. This behavior is due the fact that in this region, 
the maximum of the LLF does not lie in the vicinity of the true parameter and the 
estimator is biased. 

The cost function (6.5) is plotted in Fig. 6.2 for different observation times. The 
plots show that as the observation time increases, the cost function becomes larger. 
Figure 6.2 also show that each of the cost functions has multiple maxima. Hence, 
they verify that, to avoid getting trapped in local maxima, adopting a grid search op- 
timization approach is necessary. As the search is performed over the whole interval 
(0, 1) cycle, the algorithm initialization issues do not arise. 

It is also verified that A (0o) in (6.3) can safely be dropped from the cost function. 
The ML phase estimates are found using (6.3) as the cost function. No meaning- 
ful difference is observed between the obtained RMS errors and the ones obtained 
using (6.5). 
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Fig. 6.1 RMS error for maximum-likelihood estimator (MLE) 
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6.7.2 Pulse Delay Estimator 



Assuming the velocity data is available, the pulse delay estimator is numerically 
evaluated. The spacecraft velocities are assumed to be vi = 3 km/s and V2 = 9 km/s; 
and the initial phase observed at the first detector is (pi = 0.2 cycle. The relative 
distance between the spacecraft is assumed to be 180 km. Hence, the pulse delay to 
be estimated is t^ = 0.6 ms. 

Figure 6.3 shows the RMS of the proposed pulse delay estimator, obtained 
through simulation, against the CRLB, calculated as in (4.44). The plots show that 
as time goes on, the estimator asymptotically attains the CRLB. The cost functions 
(6.5) are optimized using a grid search approach in the domain of (0, 1) cycle. Simi- 
lar to NLS and CC estimators, as the observation time is reduced, there is a threshold 
point at which the variance of the estimation error starts to deviate from the theoreti- 
cal value. This is due to the fact that the estimator becomes biased if the observation 
time is not long enough. The reason is that as the observation time goes below this 
threshold, the number of detected photons is not high enough. Therefore, the real- 
ization of TOAs does not effectively represent the rate function of received photons. 
Hence, it results in a distorted cost function whose maximum does not lie in the 
vicinity of the true phase delay. 
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Fig. 6.3 RMS error for ML-based delay estimator 



The performance of the proposed estimators is also studied when the absolute 
velocity data is not perfectly known. The errors are assumed to be Avj = 1 ,000 and 
Av2 = 1,200 m/s. As Fig. 6.4 shows the RMS errors increase compared to the case 
where the velocity data is perfectly known. 



6.7.3 Computational Cost 



Figure 5.1 shows that asymptotic performance of the CC and NLS estimators is 
very close to the CRLB. This motivates the numerical implementation study of each 
estimator. Hence, the CPU time used by MATLAB to perform the calculations for 
one Monte Carlo realization is investigated. The utilized processor is an Intel 2.4 
GHz dual core. 

Figure 6.5 shows that the computational cost for linear interpolation used to cal- 
culate the CC estimator is almost independent of the observation time. This is to be 
expected, as the interpolation is performed at one point just to find a finer maximizer 
for the cross-correlation function. It also shows that the linear interpolation used for 
calculation of the MLE and the epoch folding have almost the same computational 
complexity, and as expected, the amount of calculations almost linearly grows as 
the observation time becomes longer. 



6.7 Numerical Examples 
„-1 
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Fig. 6.4 RMS error for ML-based delay estimator in presence of spacecraft velocity errors 
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Fig. 6.5 CPU time used for epoch folding and interpolations 
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Fig. 6.6 CPU time used for construction of the cost functions 



From Fig. 6.6, it can be seen that the amount of calculations to construct the 
CC and NLS cost functions does not change with the observation time. However, 
the computational cost to form the ML cost function almost linearly grows as the 
observation time becomes longer. It also shows that using Fourier domain approach 
to calculate the CC cost function significantly reduces the amount of calculations. 

The total CPU time used to find the CC estimator is the sum of the epoch fold- 
ing time, the time used for construction of the CC cost function, and the parabolic 
interpolation time. The total CPU time to find the NLS estimator is the sum of the 
epoch folding time and the cost function construction time for all grid points. The 
total CPU time for the MLE is the sum of the construction time and the total in- 
terpolation time for all grids. The reason is that for each grid, an interpolation is 
needed for the evaluation of the cost function at the measured TOAs. These plots 
are all shown in Fig. 6.7. 

As the plots show, the ML CPU time is bigger than the one for the NLS algorithm, 
and it grows significantly faster with the observation time. The ML calculational 
cost rate with respect to the observation time is almost 15 times bigger than the 
NLS rate in this simulation scenario. Furthermore, the NLS CPU time is about 10 
times bigger than the CPU time used to calculate the CC estimator. 



6.8 Summary 
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Fig. 6.7 Total CPU time used for calculation of the estimators 



6.8 Summary 



Using the pdf of the photon TOAs we construct an ML criterion, and formulate 
an asymptotically efficient estimator of the initial phase and the Doppler frequency 
for each detector. When the spacecraft velocity is known, it becomes only a phase 
estimator. Although the MLE asymptotically achieves the CRLB, we show that its 
computational cost is significantly higher than the epoch folding-based estimators. 
Under certain conditions, the pulse delay estimator can be used even when the space- 
craft absolute velocities are not perfectly known. We also verify the analytical results 
on the MLE's performance and the computational costs via numerical examples. 



Chapter 7 

Recursive Estimation 



7.1 Introduction 

In this chapter, a recursive algorithm is formulated which can be used to find the 
relative navigation solution between the two spacecraft. The navigation system is 
equipped with IMUs which provide the spacecraft acceleration data. The dynamics 
of relative position between the two spacecraft and a model of the IMU accuracy are 
utilized for developing the navigation algorithm. The measurements, which are ob- 
tained by time tagging the photons, are modeled as a linear function of the projected 
relative position onto the unit direction pointing to the pulsar plus the measurement 
noise. The measurement noise variance is selected based on how well the pulse de- 
lay is estimated. Then, by applying a Kalman filter, the relative position, the relative 
velocity, the relative bias between accelerometers, and the differential time between 
clocks are estimated, and the steady state estimation error covariance is obtained. 
The effect of different system parameters on the achievable accuracy of relative po- 
sition estimation is investigated. In particular, the effect of different values of IMU 
uncertainty, measurement noise variance, and number of pulsars used for estimation 
are considered. 

We formulate the system dynamics for the navigation problem in Sect. 7.2. 
The measurement equation is presented in Sect. 7.3. The Kalman filter is given in 
Sect. 7.4. We investigate the observability of the system dynamics in Sect. 7.5. In 
Sect. 7.6, we briefly explain how the proposed approach is applicable to the absolute 
navigation problem. We offer a geometric approach for estimation of the spacecraft 
absolute velocities in Sect. 7.7. Section 7.8 provides different numerical simulation 
scenarios of the navigation system. 



7.2 System Dynamics 

X-ray pulsars are being observed over the sky map. All the measurements are per- 
formed in an inertial system whose origin is the SSB. The directional location of the 
pulsars is known. We denote this directional unit vector as H , where i identifies 
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the pulsar and j denotes the spacecraft. It is assumed that the space vehicles are 
close enough so that H^ = H^ i] = ffW, 

Let Ax 6 R 3 be the relative position, Av £ R 3 be the relative velocity, and Aa £ R 3 
be the relative acceleration between the two spacecraft, and assume they follow the 
following dynamics 

AX=AAX + BAa + Gw v , (7.1) 

where AX = (Ax T Av t ) , the matrices are 

°roL> •-(£)■ -(&) ™ 

and w v is an independent zero mean white noise process with a known power spec- 
tral density (PSD), denoted by W v 

E[w v (t)w v {t;) T }=W v 8(t-T). (7.3) 

Each spacecraft is equipped with an inertial measurement unit which measures 
its acceleration, ay 1 , as 

o£u = « ( ' ) +w« +*« J »'=1.2 (7.4) 

where w' ! ' is a white noise process with a known PSD, and ba is the accelerometer's 
bias. Let Aauyiu be 

A ^ (!) ( 2 ) /"7 «\ 

AfllMU — fl IMU _fl lMU- ( - / '- ) ^ 



Therefore, from (7.4) 



where 



AaiMU = Aa + b a + w a , (7.6) 



Aa = a W -a {2) (7.7) 



and 



b a ^b { a ] -b { a ] (7.8) 

is the relative bias. Because it slowly varies, and to keep the Kalman filter open, it 
is assumed to be a Brownian motion process, as 

b a = w b , (7.9) 

where w a and w\, are zero mean independent white noise processes with known 
PSDs, denoted by W a and W h 

E[w a (t)w a (T) T ]=W !l 8(t-T) (7.10) 

E[w b (t)w b (zf]=W b S(t-T). (7.11) 
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The relative distance is determined by integration of Aoimu> 

AXimu = A AX mv + BAa mu , (7.12) 

where AXimu = (A-%rn ^ v jmu) • Defining the error as, 

X e ^AX mu -AX (7.13) 

from (7.1), (7.6), and (7.12), the dynamics of X e are given by 

X c = AX C + Bb a + Bw a + Gw v . (7. 14) 

The differential time is common among all X-ray detectors. To model its slowly 
varying dynamics, and to keep the Kalman filter open (i.e., the Kalman gain does 
not converge to zero) [58], it is assumed to be a Brownian motion process, 

te=w e , (7.15) 

where w e is an independent zero mean white noise process with a known PSD, 
named W e 

E[w e {t)w e {T) T ]=W e S(t~T). (7.16) 

Collecting the unknown vectors X e , b a , and t e into the vector X as, 

X=lb a ) (7.17) 



its dynamics are given by 



where 



X = FX + w, (7.18) 

I A B 0\ 

F= 3X 6 3x3 (7.19) 

\0lx6 lx3 0/ 



and 



A 
W = 



fw v \ 

Wb 

W a 

\We/ 



(7.20) 



is an independent white noise process where 

£[w(r)w(T) T ] =<bag{W v ,W b ,Wa,W e )8(t-T), (7.21) 
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To obtain the discrete model, the system (7.18) is sampled with the sampling 
time, T s . The discrete model is as follows [59]. 



where 



and 



X{k+l) = 0X{k) + w d (k) 
= txp(FT s ) 

w d (k)= / exp(Fr)w(T)dT. 
Jo 



(7.22) 
(7.23) 
(7.24) 



Note that F 3 = 0. Therefore, 



® = I+T,F 



which equals 



's r 2 



1t2j 







tfl 0\ 



// 3 x3 T S I 

/ 3 x3 TJ 

/ 3x3 

V 



1/ 

As w{t) is zero mean, from (7.24), w d (k) is zero mean as well 

E[w d (k)}=0 



(7.25) 



(7.26) 



and 

E[co d (k)(o d (j) T } = Q8 k] 

where Q, which is the variance matrix of the process noise, is given in (7.29). 

2 mi 

T.W, 



(7.27) 
(7.28) 



Q 



T ±w a 



T 4 






V 







^w b 





T s W h 
T s W e J 



(7.29) 



7.3 Measurements 



The pulse delay estimate for the rth pulsar is denoted by i^'. It is related to the 
relative position via 



c1®(k)=C$+Cte(k) + Tl®(k) 

= H^Ax(k)+ct e (k) + r] d i \k), 



(7.30) 



7.3 Measurements 
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where H^'> is the directional unit vector pointing to the r'th pulsar, and r\j (k) is an 
independent zero-mean white noise sequence. The pulse delay estimates are avail- 
able every T m seconds. Hence, the autocorrelation function of r\^ (k) is 



E[r\f(h)T\fV) ]=R®8 kj , /?(0 = c*var[&] (0 



(7.31) 



where var[fd]^ is the variance of the r'th pulse delay estimator, obtained every T,„ 
seconds. 

Let's define the measurement with respect to the r'th pulsar as 



yU(k)±ctH\k). 
Using (7.30), the measurements have the following discrete model 
yW (k) = H®Ax(k) + ct e (k) + r]f (k). 



(7.32) 



(7.33) 



Employing N different pulsars, N measurements are available, which can be col- 
lected into the following matrix format, 



Y(k)=CAX(k) + ctt e (k) + ri d (k), 



(7.34) 



where 



,(2) 



\yWj 



C± 



[hW 1x3 \ 
7/( 2 ) lx3 



\HW 0l J 



(7.35) 



and 



* = M" if 



»f>) 



(7.36) 
The measurement noise, r\d{k), is an independent white zero-mean process, and 

E[n d (k)T} d (j) T }=R8 kj , (7.37) 

where 

R = diag (r® } . i=\ N (7.38) 



is the measurement noise variance matrix. 

Recall that from (7.17), AX = AXtmu — X e . Therefore, (7.34) can be restated as 



Z(k)^-(Y(k)-CAX mv (k)) 
= CX e (k)-ctt e (k) + ri(k), 



(7.39) 
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where 7] (k) — —f]d{k), which is clearly an independent zero mean white noise with 
the same power spectral density as r) l i(k). Using (7.17), the measurement (7.39) can 
be expressed as a function of X(k) 

Z{k)=HX{k) + r\{k) (7.40) 

where, 

H = (C Nx3 -ct) . (7.41) 

To estimate the relative position, the relative velocity, the relative biases, and the 
differential time between clocks, the dynamic system in (7.22), and the measure- 
ment equation in (7.40) will be employed in a Kalman filter. 



7.4 Discrete Time Estimation Process 

Let X k be the a priori system state estimate at time k, and P k be the a priori 
estimation error covariance, which are based on all measurements up to Z k _ i . Sim- 
ilarly, let X k and P k denote the a posteriori estimate and covariance matrix after 
the measurement Z k has been processed. Assuming that measurements start at k = 0, 
the Kalman filter equations in discrete time are as follows [58,59]. 
The filter's initial values are set as 

X H = E[Xo] (7.42a) 

P H =E[(X -x{ ) - ) )(Xo-x{ ) - ) ) T } (7.42b) 

and the state estimates are updated as 

K k = P k H H T [HP { k ~ ] H T + R}- 1 (7.43a) 

X i +) = X t ] + K klZk - HX { - ] ] (7.43b) 

P k i+) = [I- K k H]P k ( - ] [I - K k H] T + K k RKj (7.43c) 

X i+l = 0X k +) ( 7 - 43d ) 

p k+l = ®P { k +) $ T + Q- (7.43e) 

Note that the measurements (7.40) are available at a lower rate than the posi- 
tion estimate updates. Hence, when no measurement is available, the a priori and a 
posteriori estimates are equal. In other words, (7.43b) must be replaced by 

x { k +) = x t ] ( 7 - 44 ) 

in the meantime between the pulsar measurements (7.40) become available. 



7.5 Discussion 

7.5 Discussion 
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For the Kalman filter to estimate all the unknown states, the dynamic system (7.22) 
and (7.40) must be observable. Hence, the observability matrix is investigated. 
It equals 



HF 
HF 2 J 




r o o 

oro 
o o r 


-ct 





(7.45) 



where 



n 



Afx3 



H (2) 



\ H W) 



(7.46) 



Note that & is an 3N x 10 matrix. If there are at least four different pulsar mea- 
surements (N > 4) where no two are along the same directional vector from the 
spacecraft, then & will be full column-rank (i.e., rank 10). Hence the continuous 
system (7.18) is observable using the measurement (7.40). To conclude about the 
observability of discrete dynamics (7.22) and (7.40), the following theorem is em- 
ployed [60]. 

Theorem 7.1. Suppose the continuous system (7.18) and (7.40) is observable. A 
necessary and sufficient for its discretized equations (7.22) and (7.40), with sam- 
pling period T s , to be observable is that |3m[A,- — Xj\ \ ^ 2nm/T s for m — 1,2,..., 
whenever !SHe[A; — Xj\ = 0, where Xj is an eigenvalue ofF. 

Because all eigenvalues of the matrix F are zero, then 



\3m[Xi-Xj}\ = 



+ 



2nm 



(7.47) 



Therefore, it can be concluded that the discretized system is fully observable too, if 
6 is full column-rank. 

As a result, by choosing at least four different pulsars with different direction 
vectors, all the states can be estimated using the Kalman filter. 
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7.6 Absolute Navigation 

The proposed navigation approach is applicable for absolute navigation as well. It 
is clear that if the location of one of the spacecraft is known, the relative naviga- 
tion problem boils down to an absolute navigation problem. As the location of the 
SSB is known, placing one of the spacecraft at the SSB, using the same problem 
formulation, the absolute navigation solution can be obtained. 



7.7 A Geometric Approach for Estimation of Absolute Velocities 

In this section, we propose a different formulation of the navigation problem, using 
geometric dispersion of the X-ray pulsars over the sky map. 

For simplicity, a two-dimensional scenario where four pulsars are utilized is con- 
sidered (see Fig. 7.1). The ith pulsar's estimated time delay is composed of the true 
pulse delay and the differential time between clocks. From (5.1), and using Taylor 
series, it equals 



i (,) -t ( 






{l+V J cos9,/c)f i 



ti' 



,(fl 



(0 

Vj cos Qi 



(7.48) 




Fig. 7.1 Using pulsar geometry to estimate absolute velocities 
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where A0 ^ = 0i ' — 02 ' , Vj is the absolute velocity vector of the jth spacecraft, 0,- 
is the angle between Vj and H®, and m ! ' is the uncertainty modeling for the phase 
estimation error. Hence, 

-V = fi ' +t e + -fc-Vj cos 0,- + tj W , (7.49) 

where 

,» = - M M(V^°^Y (7.50) 



Defining the (th pulsar measurement as jW = A0W// s , (7.49) becomes 

z C0=#> + (<+z (0^^ + T ,(0. C7.51) 

c 

Because fl = H^'Ax/c, the new set of measurements is 

z (,) = +f f + z w -^ ! - + T7 (0 . (7.52) 

c c 

Furthermore, the angular relations between directional vectors form a set of 
pseudo-measurements 

Of = 61-62+71 (7.53a) 

91 = 2 + 03 + 72 (7.53b) 

0P-04-03 + 73, (7.53c) 

where Jt is the uncertainty. Using (7.52) and (7.53) as measurements, a new Kalman 
filter can be formulated for estimation of Ax, V,-, and 0,. 



7.8 Numerical Examples 

To simulate the three-dimensional navigation algorithm, eight pulsars are selected 
(see Fig. 7.2 for pulsar geometry). The pulsar profiles are shown in Fig. 4.2. The 
pulsar galactic coordinates, and their photon flux values are given in Table 7. 1 [19]. 

The detector area is assumed to be 10 4 cm 2 . The effective pulsar rate, A s , is 
calculated according to the pulsar photon flux and the detector area. The effective 
background rate is chosen as Ab = [As/10] . The corresponding values are given in 
Table 7.2. 

In three different scenarios, the measurements are incorporated every 100 s, 
and every 500 s. The NLS-based pulse delay estimation approach is utilized, and 
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Fig. 7.2 Pulsar geometry on the sky map 

Table 7.1 Employed pulsars 



Pulsar 


Period (s) 


Galactic 
longitude (o) 


Galactic 
latitude (o) 


Flux (2-10 keV) 
(ph/cm 2 /s) 


B0531+21 


0.0335 


184.56 


-5.78 


1.54E + 00 


B0540 - 69 


0.0504 


279.72 


-31.52 


5.15E-03 


B0833 - 45 


0.0893 


263.55 


-2.79 


1.59E-03 


B1509-58 


0.1502 


320.32 


-1.16 


1.62E-02 


B1821-24 


0.0031 


7.80 


-5.58 


1.93E-04 


B1937 + 21 


0.0016 


57.51 


-0.29 


4.99E - 05 


B1055-52 


0.1971 


164.50 


-52.45 


1.64E-06 


J0437 -47 


0.0057 


253.39 


-41.96 


6.65E - 05 


Table 7.2 Employed pulsars 


(detector area: 


10 4 cm 2 ) 




Pulsar 


A b 
(ph/s) 


A s 
(ph/s) 


<T„, (m) 
(_T obs = 100s) 


a m (m) 

(r obs =500s) 


B0531+21 


1.54E + 03 


1.54E + 04 


5.64E + 2 


2.52E + 2 


B0540 - 69 


5.15E + 00 


5.15E + 01 


3.46E + 4 


1.54E + 4 


B0833 - 45 


1.59E + 00 


1.59E + 01 


4.67E + 4 


2.09E + 4 


B1509-58 


1.62E + 01 


1.62E + 02 


8.62E + 4 


3.86E + 4 


B1821-24 


IE + 00 


1.93E + 00 


4.56E + 3 


2.04E + 3 


B1937 + 21 


IE + 00 


4.99E-01 


4.53E + 3 


2.02E + 3 


B1055-52 


IE + 00 


1.64E-02 


1.67E + 8 


7.49E + 7 


J0437 -47 


IE + 00 


6.65E - 03 


1.11E + 5 


4.99E + 4 



(0 



Om = V/?W values corresponding to each scenario are given in Table 7.2. The 
spacecraft absolute velocity data are not perfectly known. The velocity errors are 
assumed to be Avi = 1000 and Av2 = 1200 m/s. Note that because of the effect of 
absolute velocity errors on the quality of measurements, the measurement noise 
variance values are bigger than their corresponding analytical values when the 
velocity data is perfectly known. 
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To initiate the Kalman filter, the a priori covariance matrix is chosen as 

P = dmg{P 0x ,Po v ,Po ba ,Po te ), (7.54) 

where 

V^) A - = diag(10 6 ,10 6 ,10 6 ) (m) (7.55a) 

V^ v = diag(10 4 ,10 3 ,10 2 ) (m/s) (7.55b) 

^ fl =diag(10 a5 ,10 - 5 ,10 a5 ) (m/s 2 ) (7.55c) 

/^ = 10 3 (s). (7.55d) 

The a priori state estimate at k = is Xq ' = 0. Note that although the initial 
uncertainties given in (7.55) change the transient response of the Kalman filter, they 
do not affect its steady state performance. 

The IMU process noise PSDs used by the Kalman filter are given in (7.56). Note 
that to simulate the dynamics (7.22), it is assumed that Wj, = and W e = 0, since 
they are assumed to be constant states. Although, they follow slowly time varying 
dynamics in practice. Hence, to take this fact into account, and to keep the Kalman 
filter open, they are modeled as Brownian motion processes in the estimation stage. 

(7.56a) 
(7.56b) 
(7.56c) 
(7.56d) 

The units in (7.56) are obtained from (7.29), noting that the diagonal elements 
have units m 2 , (m/s) 2 , (m/s 2 ) 2 , and s 2 , respectively. 

In (7.56), the accelerometer bias uncertainty, Wj,, is chosen small based on the 
fact that the available IMUs are very accurate. These values are easily achievable 
by available commercial IMUs such as LN200 [61]. Choosing other W\, values 
only affects achievable estimation accuracies. The accelerometer uncertainty, W a , 
represents the applied forces to the spacecraft from the celestial sources (ex. solar 
radiation pressure). Hence, it is chosen small as well. 

To find out how accurate the relative distance between the spacecraft can be esti- 
mated, the state estimation error, and the standard deviation (STD) of the a posteriori 
estimation error are investigated. In other words, the square roots of the diagonal el- 
ements of P^ are considered. Note that using (7.17), the standard deviation of X e 
equals to the standard deviation of AX, which represents the relative position and 
relative velocity vectors between the spacecraft. All results are obtained through 
Monte Carlo simulation over 3,000 independent realization of the stochastic signals. 
Different sampling times, T s , are chosen for discretizing the continuous dynamics. 
As expected, T s does not change either the error covariance in the steady state or the 
time that the steady state is reached; however, the number of iterations to get to the 
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same values of the error covariance varies. The initial state of the dynamic system 
(7.22) is assumed to be a Gaussian random variable 



Xb~^(0,/b), 



(7.57) 



where Pq is given in (7.54). 

In all scenarios the Monte Carlo state estimation error, obtained by averaging 
over 3,000 realizations, is observed to be zero. Hence, its plots are omitted. The 
Monte Carlo a posteriori standard deviation, along with the analytical one obtained 
from (7.43c) are plotted, and it is shown that they satisfactorily match. The analyt- 
ical values that are computed in the Kalman filter are based on modeling some of 
the variables as Brownian motion, as opposed to the simulation model where the 
variables are assumed to be slowly varying (in some cases constants). As a result, 
the simulation STD plots may not perfectly match the analytical values, but they are 
very close. 

The estimation error obtained from one [example] realization is plotted as well. 
As expected, we can see that at any time epoch, most of the Mote Carlo estimation 
error points are within the 1-CT envelope. 

In the first scenario, eight pulsars listed in Table 7.1 are selected, and the mea- 
surements are obtained every 100 s. The results are plotted in Figs. 7.3-7.6. The 
obtained estimation accuracies are given in Tables 7.3, 7.4, and 7.5. 
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Fig. 7.6 Clock differential time estimation (N = 8, r obs = 100 s) 
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Table 7.3 Position 
estimation results 



Scenario 



fPAm) 



fP y (m) 



n = %, r obs = ioos 

N = 8, r obs = 500 s 
A' = 4,7 obs = 100s 



/^(m) 



1.93E + 3 6.74E + 3 3.87E + 4 
1.16E + 3 4.97E + 3 3.21E + 4 
2.16E + 4 4.12E + 4 6.83E + 4 



Table 7.4 Velocity 


estimation results 






Scenario 


\/P^(.m/s) 


V^Cm/s) 


y/\ (m/s) 


iV=8,r obs = 100s 


7.52 


21.87 


92.96 


N = 8, r obs = 500 s 


2.09 


7.91 


47.55 


iv=4,r obs = ioos 


15.21 


114.78 


183.08 



Table 7.5 IMU Bias and clock differential time estimation results 




Scenario 


sJP^ (m/s 2 ) 


\AF(m/s 2 ) 


y^pCm/s 2 ) 


VE (sec) 


/V = 8, r obs = 100s 


1.31E-2 


3.74E-2 


1.34E-1 


1.32E-5 


/V = 8,r obs = 500s 


1.59E-3 


5.87E-3 


3.45E-2 


1.14E-5 


N = 4, r obs = 100s 


2.05E-2 


1.80E-1 


2.15E-1 


8.01E-5 
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In the second scenario, employing the same pulsars, the measurements are 
provided every 500 s. Hence, the measurement noise STD values decrease. The 
corresponding values are given in Table 7.2. As Figs. 7.7-7.10 show, using more 
accurate measurements, results in obtaining more accurate estimates. The estima- 
tion results are given in Tables 7.3, 7.4, and 7.5. 

In the third scenario, to investigate the effect of the pulsar geometry on the esti- 
mation accuracy, just the first four pulsars from Table 7.1 are chosen for navigation. 
The measurements are updated every 100 s. The plots are shown in Figs. 7.1 1-7.14, 
and the resulted accuracies are given in Table 7.3. As expected, compared to the 
case where eight pulsars were employed, the estimation error standard deviations 
have increased (Tables 7.3, 7.4, and 7.5). 

Recall from (7.33) that the relative position, Ax, and the differential time between 
clocks, t e , show up directly in the measurements. Hence, we expect the Kalman 
filter to be able to estimate them immediately, and it takes longer for the relative 
bias and velocity estimates to converge to the steady state values. From the a priori 
initial standard deviations, given in (7.55), and using the estimation plots shown in 
Figs. 7.3-7.14, this can be verified where the a posteriori position and differential 
time estimation standard deviations dramatically drop when the first measurements 
is incorporated at t = 0. 

Another interesting point regarding Tables 7.3, 7.4, and 7.5 is that the 
x-coordinate estimates are more accurate than the y-coordinate estimates. Consider- 
ing Fig. 7.2, we expect this phenomenon. As the sky map shows, the employed 
pulsars are geometrically more distributed along the jc-axis than the y-axis. 
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Fig. 7.7 Relative position estimation (N = 8, r obs = 500 s) 
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Fig. 7.8 Relative velocity estimation (N = 8, T obs = 500 s) 



Furthermore, Tables 7.3, 7.4, and 7.5 show that the z-coordinate estimates are 
less accurate than the x- and y-coordinate estimates. This implies that the pulsars 
geometric distribution along the z-axis is worse than the other axes. 

In summary, we investigated effectiveness of the proposed algorithm through 
different simulation scenarios. Achievable estimation accuracies depend on different 
elements of the navigation system. The main factors are the pulse delay estimation 
accuracies, the IMU accuracies, and the pulsars geometric distribution on the sky 
map. 

We have also performed a parametric study to investigate the range of achievable 
accuracies as a function of the quality of measurements. All simulation parameters 
are the same as the scenario where N = 8 and r o b s = 500 s. The only difference 
is the steady state STD of estimation errors are obtained for kR<S m in the range of 
0.1 < Icr < 10, where a m values are the measurement noise standard deviations for 
?obs = 500 s, given in Table 7.2. The results are plotted in Figs. 7.15-7.18. From 
these plots, the range of estimation accuracies are given in Table 7.6. 

These ranges show that employing bright pulsars, position estimation accuracies 
in the order of a couple of 100 m are achievable for position estimation. The velocity 
estimation accuracies in the order of a few meters per second can be obtained. The 
attainable bias estimation accuracy is less than 1 mm/s 2 . Furthermore, the differen- 
tial time between clocks can be estimated in the order of a few micro seconds. 

The effect of the IMU uncertainty on relative navigation solution is also inves- 
tigated. All parameters are the same as the scenario where N = 8 and r o b s = 500 s. 
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Fig. 7.9 Relative bias estimation (TV = 8, T obs = 500 s) 



4 


X 


1<T 




















1 ' 














3 












Analytic STD 

- Monte Carlo STD 


- 


2 










- 








- 


- 




1 

IT 
J2, 


'""'""-. ...^ 


-■'■''""'--. 


-1 










- 




'■""' ■ * 




-2 










- 


-3 

4 








1 1 


1 1 1 1 



10 15 20 25 30 35 

Time (min) 



40 



Fig. 7.10 Clock differential time estimation (N = 8, r obs = 500 s) 
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Fig. 7.11 Relative position estimation (JV = 4, 7/ obs = 100 s) 
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18 



18 



** — ^^^ \ y 'm- ■ ' " ' 

N. • "*■ ' i i 



18 



7.8 Numerical Examples 



105 



S 



-5 



ri" 



12 



Time (min) 
Fig. 7.13 Relative bias estimation (TV = 4, r obs = 100 s) 
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Fig. 7.15 The STD of position estimation error 




Fig. 7.16 STD of velocity estimation error 
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Fig. 7.17 The STD of bias estimation error 




Fig. 7.18 The STD of differential time estimation error 
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Table 7.6 Estimation results when R or Q matrices change 



R varies 



Min 



Max 



<2 varies 



Min 



Max 



/K' 0*0 

/P~(m/s) 
/^(m/s) 
/PT(nVs) 

/iJF (m/s 2 ) 
/p|"> (m/s 2 ) 
/p^ (m/s 2 ) 



142.42 
516.83 
3051.2 

0.2489 
0.5076 
1.9286 

3.3694E-4 
4.1967E-4 
6.9230E-4 

1.1508E-6 



8438.4 
2946.8 
169950 

4.6571 

13.6063 

64.1779 

0.0012 
0.0032 
0.0137 

5.7419E-5 



852.27 
2947.5 
17137.0 

0.4737 
1.3696 

7.0152 

0.12515 
0.3271 
0.0016 
5.7494E-6 



1423.6 
5157.0 
30406.0 

2.4889 
5.0678 
19.1681 

0.0034 

0.0042 

0.0069 

1.1473E-5 
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Fig. 7.19 The STD of position estimation error 



The only parameters that are changed are the IMU uncertainties. Let o p be any of the 
IMU uncertainties given in (7.56d). Then, the achievable estimation accuracies are 
examined for kgO p , where 0.1 <kg< 10. The results are plotted in Figs. 7. 19-7.22. 
The minimum and maximum accuracies are also given in Table 7.6. 
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Fig. 7.20 The STD of velocity estimation error 




Fig. 7.21 The STD of bias estimation error 
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Fig. 7.22 The STD of differential time estimation error 



7.9 Summary 



In this chapter, we present a recursive navigation algorithm. We utilize the pulse 
delay estimates to construct the measurements for a Kalman filter. We suggest to 
employ IMUs on each spacecraft to provide the acceleration data. Using the Kalman 
filter, we offer a recursive algorithm to estimate the relative position, the relative ve- 
locity, the relative accelerometer biases, and the differential time between clocks. 
Through different simulation scenarios, applicability of the propose navigation al- 
gorithm is also verified. 



Chapter 8 
Epilog 



This book has proposed a new approach for navigation of spacecraft in space 
employing X-ray pulsar measurements. The presented approach is applicable for 
both absolute and relative navigation problems. Pulsars emitting in the X-ray band 
were chosen because of their stable period and their geometric distribution in the 
sky map. The main advantage of using X-ray pulsars for navigation is that rela- 
tively small size detectors can be used for detection of the X-ray photons on board 
a spacecraft. This facilitates the spacecraft design procedure. 

The developed navigation technique is based on utilizing X-ray detectors on each 
spacecraft and locking the detectors on the same pulsar. Hence, the vehicle farther 
from the pulsar detects a signal whose intensity is the time delayed version of the 
one detected by the closer vehicle. The distance between the space vehicles is pro- 
portional to the time delay. The proposed approach is to periodically estimate the 
pulse delay, and then use these estimates as measurements in a recursive algorithm 
to find the navigation solution. To analyze the system, mathematical models were 
developed for the X-ray pulsar signals, and the Cramer-Rao lower bound for es- 
timation of the pulse delay was presented. Two different strategies for estimation 
of the pulse delay were suggested. One strategy was to employ epoch folding. 
The procedure, which is used to recover photon intensities, was mathematically 
studied and characterized. Two different estimators based on epoch folding were 
formulated and their asymptotic performance was analyzed. One estimator uses the 
cross-correlation function between the empirical rate function and the known pul- 
sar intensity function. The other estimator is obtained through solving a nonlinear 
least squares problem. It was shown that these estimators are consistent, but not 
asymptotically efficient. The second strategy was to directly utilize the measured 
photon time of arrivals. Using this strategy and based on a maximum likelihood cri- 
terion, a pulse delay estimator was formulated, and it was shown that the estimator 
is asymptotically efficient. It was shown that the cross-correlation-based estimator 
is computationally more efficient than the other two estimators. 

Space vehicles are equipped with inertial measurement units to provide the accel- 
eration information, which are converted to velocity data and utilized by the pulse 
delay estimator. The pulse delay estimates, in turn, are taken in as measurements 
by the Kalman filter for recursive estimation. The measurement noise variance 
is selected based on accuracy of the pulse delay estimates. Models of spacecraft 
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dynamics and inertial measurement units are employed by the proposed algorithm. 
It was shown that the relative inertial measurement unit biases and the differen- 
tial time between detectors' clocks can be estimated as well as the relative position 
and velocity between the spacecraft. The three-dimensional relative navigation so- 
lution can be obtained by taking measurements on four or more different X-ray 
pulsars. Depending on the number of pulsars, their characteristics, their geometric 
distribution in the sky, and inertial measurement unit uncertainties, the achievable 
estimation accuracies were examined. 

To enhance the position estimation accuracy, improving models of the spacecraft 
motion and IMU dynamics is necessary for future work. Because some models are 
now nonlinear, employing filtering techniques such as the extended Kalman filter 
and the unscented Kalman filter may be appropriate. Modeling clock errors and 
studying their effect on the navigation solution are necessary to obtain more accurate 
results. Furthermore, to estimate the system biases with more precision, smoothing 
algorithms may be utilized. Developing signal processing techniques for estimation 
of the pulse phase when spacecraft velocities are varying rapidly is another necessity 
for improving the navigation algorithm. Additionally, it is of interest to investigate 
the navigation problem in more detail in situations where Doppler frequencies and 
pulse phases are simultaneously estimated. Proposing pulse delay estimators for 
these situations and analysis of their performance are important areas of future re- 
search. Addressing the cycle ambiguity problem is another topic to work on. Also, 
taking advantage of pulsar geometry for new formulations of the navigation prob- 
lem, is an interesting field for enhancing the pulsar-based navigation solution. 
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